bge-barcoding · SchistoDan · Aug 21, 2026 · Aug 19, 2026 · Aug 19, 2026 · Aug 19, 2026
diff --git a/CHANGELOG b/CHANGELOG
@@ -1,5 +1,33 @@
 # Changelog
 
+## Unreleased
+
+- Taxonomic validation: `tv_local_blast.py` now batches queries into chunks of 50
+  sequences per `blastn` invocation instead of running one process per sequence,
+  so a large reference database is loaded once per chunk rather than once per
+  barcode. Per-sequence TSV outputs, the summary CSV layout and resume behaviour
+  are unchanged. A failed chunk retries its sequences individually, and
+  sequences that still fail now cause a non-zero exit (with no summary CSV)
+  rather than being reported as no-match; `--allow-partial-failures` restores the
+  previous leniency.
+- Fixed the summary CSV dropping all hits for sequences whose FASTA header
+  carries a description, where BLAST's `qseqid` (first token only) did not match
+  the sanitized full header used as the CSV key.
+- Removed the unused `blast_options` entry from the Snakefile's
+  `taxonomic_validation` fallback, and corrected the documented BLAST hit limit
+  (100, not 500).
+- Structural validation: fixed a temp-file leak in `structural_validation.py`,
+  where the per-sequence nhmmer query and tabular files were created with
+  `delete=False` and never removed, leaving two files per sequence in `TMPDIR`
+  (48 per sample, since each sample is validated across 6 parameter
+  combinations and 4 input FASTAs).
+- Structural validation: added `--threads`, passed from the rule, setting
+  nhmmer's `--cpu` per invocation. Defaults to 1, so standalone behaviour is
+  unchanged.
+- Right-sized the `structural_validation` rule resources from 32 GB / 1 thread /
+  himem to 4 GB / 8 threads / medium. Peak memory scales with sequence count
+  rather than read depth, measured at roughly 110 MB per 1000 samples.
+
 ## v3.0.4 (2026-08-13)
 
 Initial PyPI release.

diff --git a/README.md b/README.md
@@ -137,7 +137,7 @@ bash run_local.sh
    - Cleaned consensus generation and metrics aggregation ([`05_consensus_generator.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/05_consensus_generator.py))
 6. **Barcode validation and selection** (see [Validation process](#validation-process)):
    - **Structural validation** - HMM-based barcode extraction, reading frame analysis, stop codon detection and quality ranking ([`structural_validation.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/structural_validation.py))
-   - **Local BLASTn search** - parallel searches of structurally validated barcodes against a local reference database ([`tv_local_blast.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/tv_local_blast.py))
+   - **Local BLASTn search** - chunked, parallel searches of structurally validated barcodes against a local reference database ([`tv_local_blast.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/tv_local_blast.py))
    - **Taxonomic validation** - hierarchical matching of BLAST results against expected taxonomy, selecting the best sequence per sample ([`tv_blast2taxonomy.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/tv_blast2taxonomy.py))
 7. **Statistics compilation** - QC, recovery, cleaning, filtering and validation metrics aggregated into CSV reports ([`compile_barcoding_stats.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/compile_barcoding_stats.py)).
 8. **Final integration** - all pipeline metrics merged into a unified output CSV ([`val_csv_merger.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/val_csv_merger.py)).
@@ -504,7 +504,7 @@ Structural validation (via `structural_validation.py`) assesses every barcode co
 ## Taxonomic validation ##
 Taxonomic validation runs in two steps, via `tv_local_blast.py` and `tv_blast2taxonomy.py`.
 
-**1. Local BLASTn search.** Parallel BLASTn searches against a local database, either built from a multi-FASTA with `makeblastdb` or supplied pre-built. The e-value threshold is hardcoded to 1e-5. Per-sequence TSV outputs (outfmt 6) hold the top 500 hits ordered by descending percent identity; the top 100 are carried into the summary CSV.
+**1. Local BLASTn search.** BLASTn searches against a local database, either built from a multi-FASTA with `makeblastdb` or supplied pre-built. The e-value threshold is hardcoded to 1e-5 and `max_target_seqs` to 100. Queries are batched into chunks of 50 sequences per `blastn` invocation, with up to `threads` chunks running concurrently, so the database is loaded once per chunk rather than once per sequence. Results are demultiplexed back into one TSV per sequence (outfmt 6), each holding up to 100 hits ordered by descending percent identity, and all 100 are carried into the summary CSV. Sequences whose TSV already exists are skipped, so an interrupted run resumes without repeating finished work.
 
 **2. Taxonomic assignment validation.** BLASTn results are checked against expected taxonomy using hierarchical matching and quality-based filtering:
 1. Parse the local BLASTn summary CSV, per-sample expected lineages, database taxonomy mappings, and structurally validated sequences.

diff --git a/beegees/config/config.yaml b/beegees/config/config.yaml
@@ -97,72 +97,81 @@ taxonomic_validation:
 # taxdump_dir: "/path/to/writable/ncbi_taxdump"
 
 
-## Resource allocation for each rule. Memory shown in Mb (mem_mb), e.g. 2048 = 2G memory. Rules have dynamic memory scaling upon retry (mem_mb * retry #).
+## Resource allocation for each rule
+#
+# mem_mb is in MB (2048 = 2 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.
+#
+# 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
+# highlight what drives each one.
+#
+# partition is only read for gene_fetch, MitoGeneExtractor, structural_validation
+# and taxonomic_validation. Every other rule takes the profile default (medium).
+# Do not delete those four keys: the Snakefile indexes them directly and will fail
+# with a KeyError.
 rules:
-  gene_fetch:                                    # Gene Fetch: reference retrieval from NCBI GenBank
-    mem_mb: 8192
-    threads: 4
-    partition: long
+  # ===== PER-SAMPLE: cost per job is independent of how many samples are in the run =====
+  gene_fetch:                                    # gene-fetch: reference retrieval from NCBI GenBank
+    mem_mb: 4096
+    partition: PLACEHOLDER                       # Runtime scales on sample number. 
   fastp_qc:                                      # fastp: adapter trimming, QC, merging/deduplication (PE merge and concat modes)
-    mem_mb: 16384
-    threads: 4
-  clean_headers_merge:                           # Clean fastq headers for MitoGeneExtractor compatibility (merge mode)
     mem_mb: 8192
-    threads: 1
+    threads: 4                                   # fastp --thread. Note gains flatten past ~8 even on deep samples
+  clean_headers_merge:                           # Clean fastq headers for MitoGeneExtractor compatibility (merge mode)
+    mem_mb: 4096
   fastq_concat:                                  # Concatenation of fastp QC'd R1+R2 fastq files (concat mode)
-    mem_mb: 8192
-    threads: 1
+    mem_mb: 2048
   quality_trim:                                  # Trim Galore: secondary quality trimming of concatenated reads (concat mode)
     mem_mb: 8192
-    threads: 4
+    threads: 4                                   # trim_galore --cores. Note it spawns ~3-4x this many processes (cutadapt + pigz), so 4 occupies ~15 CPUs
   downsample:                                    # reformat.sh: optional read downsampling (concat/merge/SE modes)
     mem_mb: 8192
-    threads: 4
-  MitoGeneExtractor:                             # MitoGeneExtractor: protein-guided barcode extraction (concat/merge/SE modes); single-threaded (Exonerate)
-    mem_mb: 32768
-    threads: 1
-    partition: himem
-  rename_and_combine_cons:                       # Standardise and combine MGE consensus FASTA headers across parameter combinations
-    mem_mb: 8192
-    threads: 4
+  MitoGeneExtractor:                             # MitoGeneExtractor: protein-guided barcode extraction (concat/merge/SE modes)
+    mem_mb: 32768                                # Per sample per r/s combination; scales with read depth.
+    partition: PLACEHOLDER                       # Submit to high memory nodes (iv available) to allow memory scaling upon retry.
+  # ===== AGGREGATE: one job for the whole run, so memory grows with sample count =====
+  rename_and_combine_cons:                       # rename_headers.py: standardise and combine MGE consensus headers across parameter combinations
+    mem_mb: 4096
+    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: 8192
-    threads: 4
+    mem_mb: 2048                                 # Flat memory requirement regardless of sample number
   human_cox1_filter:                             # 01_human_cox1_filter.py: remove human COX1 contamination
-    mem_mb: 16384
+    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
   at_content_filter:                             # 02_at_content_filter.py: remove reads with aberrant AT content
-    mem_mb: 16384
+    mem_mb: 8192                                 # As above: [threads] x [largest alignment]
     threads: 4
   statistical_outlier_filter:                    # 03_statistical_outlier_filter.py: remove reads statistically dissimilar to the consensus
-    mem_mb: 16384
+    mem_mb: 8192                                 # As above: [threads] x [largest alignment], plus per-alignment numpy distance arrays - the heaviest of the three filters
     threads: 4
   reference_filter:                              # 04_reference_filter.py: optional reference-based sequence retention or contaminant removal
     mem_mb: 4096
-    threads: 1
+    threads: 4
   consensus_generation:                          # 05_consensus_generator.py: generate cleaned consensus sequences and aggregate metrics
-    mem_mb: 8192
+    mem_mb: 8192                                 # As above: [threads] x [largest alignment].
     threads: 4
-  extract_stats_to_csv:                          # mge_stats.py: extract per-sample MGE alignment statistics to CSV
-    mem_mb: 8192
-    threads: 1
-  structural_validation:                         # structural_validation.py: HMM-based extraction, reading frame and stop-codon checks; nhmmer is single-threaded
-    mem_mb: 32768
-    threads: 1
-    partition: himem
-  taxonomic_validation:                          # tv_local_blast.py: parallel BLASTn searches against local reference database
-    mem_mb: 8192
-    threads: 8
-    partition: medium
+  extract_stats_to_csv:                          # compile_barcoding_stats.py: extract per-sample MGE alignment statistics to CSV
+    mem_mb: 4096                                 # Builds one table across all samples, so grows with sample count
+  structural_validation:                         # structural_validation.py: HMM-based extraction, reading frame and stop-codon checks
+    mem_mb: 4096
+    threads: 4                                   # CPUs given to each nhmmer call (--cpu). Processing is sub-linear, so >4 threads has diminishing returns
+    partition: PLACEHOLDER
+  taxonomic_validation:                          # tv_local_blast.py: chunked BLASTn searches against local reference database
+    mem_mb: 8192                                 # ~250 MB per concurrent (50 consensus seq) chunk, plus the summary table over all barcodes
+    threads: 8                                   # Max concurrent single-threaded BLASTn chunks (50 sequences per chunk).
+    partition: PLACEHOLDER
   blast2taxonomy:                                # tv_blast2taxonomy.py: hierarchical taxonomic matching and best-barcode selection
-    mem_mb: 16384
-    threads: 4
-  download_taxdump:                              # Download NCBI taxdump for blast2taxonomy
     mem_mb: 4096
-    threads: 1
-  multiqc_plots:                                 # Generate per-sample summary plots from pipeline metrics
-    mem_mb: 8192
-    threads: 2
+  download_taxdump:                              # Download and extract the NCBI taxdump for blast2taxonomy
+    mem_mb: 2048
+  multiqc_plots:                                 # multiqc_plots.R: per-sample summary plots from pipeline metrics
+    mem_mb: 4096                                 # R data frames spanning all samples; grows with sample count
   multiqc:                                       # MultiQC: aggregate all QC reports into a single interactive HTML report
-    mem_mb: 8192
-    threads: 1
+    mem_mb: 8192                                 # Parses every report in the run - the largest aggregate consumer at high sample counts
diff --git a/beegees/config/multiqc_config.yaml b/beegees/config/multiqc_config.yaml
@@ -6,8 +6,8 @@ title: "BeeGees Pipeline Report"
 
 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) Snakemake workflow</a>
-  developed at NHMUK. The report outlines overall barcoding success per sample, as well as read Quality Control (QC)
+  <a href="https://github.com/bge-barcoding/BeeGees">BeeGees (Barcode Gene Extraction and Evaluation from Genome Skims) pipeline</a>
+  developed at NHMUK by Dr. Dan Parsons. The report outlines overall barcoding success per sample, as well as read Quality Control (QC)
   statistics, the outcome of reference retrieval (Gene Fetch), and barcode validation results.
   This report is complemented by the <code>*_final_metrics.csv</code> and plots within <code>05_barcoding_outcome/plots/</code>.
 

diff --git a/beegees/workflow/Snakefile b/beegees/workflow/Snakefile
@@ -759,8 +759,9 @@ if downsampling_enabled:
                 echo "Placeholder output created: $(date)" >> {log.out}
                 exit 0
             fi
-
-            reformat.sh in={input.merged_reads} \\
+
+            reformat.sh -Xmx$(( {resources.mem_mb} * 80 / 100 ))m -eoom \\
+                        in={input.merged_reads} \\
                         out={output.downsampled} \\
                         samplereadstarget={params.max_reads} \\
                         sampleseed=12345 \\
@@ -1060,7 +1061,8 @@ if downsampling_enabled:
                 exit 0
             fi
 
-            reformat.sh in={input.concat_trimmed} \\
+            reformat.sh -Xmx$(( {resources.mem_mb} * 80 / 100 ))m -eoom \\
+                        in={input.concat_trimmed} \\
                         out={output.downsampled} \\
                         samplereadstarget={params.max_reads} \\
                         sampleseed=12345 \\
@@ -1186,7 +1188,6 @@ if downsampling_enabled:
         log:
             out=os.path.join(preprocessing_dir_se, "logs/downsample/{sample}.out"),
             err=os.path.join(preprocessing_dir_se, "logs/downsample/{sample}.err")
-        threads: rule_resources["downsample"]["threads"]
         resources:
             mem_mb=lambda wildcards, attempt: rule_resources["downsample"]["mem_mb"] * attempt
         retries: 3
@@ -1206,7 +1207,8 @@ if downsampling_enabled:
                 exit 0
             fi
 
-            reformat.sh in={input.se_trimmed} \\
+            reformat.sh -Xmx$(( {resources.mem_mb} * 80 / 100 ))m -eoom \\
+                        in={input.se_trimmed} \\
                         out={output.downsampled} \\
                         samplereadstarget={params.max_reads} \\
                         sampleseed=12345 \\
@@ -3064,11 +3066,11 @@ rule multiqc_plots:
             {params.taxdump_dir} \
             {params.run_mode} \
             >> {log} 2>&1
-
+        
         # Move static PNGs from mqc_dir to plots_dir
         find {params.mqc_dir} -maxdepth 1 -name "*.png" \
             -exec mv {{}} {params.plots_dir}/ \;
-
+        
         echo "Completed: $(date)" >> {log}
         """
 
@@ -3091,29 +3093,29 @@ rule multiqc:
     shell:
         """
         set -euo pipefail
-
+        
         mkdir -p {params.output_dir}
         mkdir -p $(dirname {log})
-
+        
         echo "Running MultiQC: $(date)" > {log}
         echo "Data dir:    {params.mqc_data_dir}" >> {log}
         echo "Config:      {input.mqc_config}" >> {log}
         echo "Output dir:  {params.output_dir}" >> {log}
-
+        
         multiqc {params.mqc_data_dir} \
             --config {input.mqc_config} \
             --outdir {params.output_dir} \
             --filename multiqc_report.html \
             --force \
             >> {log} 2>&1
-
+        
         # Copy report back to main output directory
         cp {params.output_dir}/multiqc_report.html \
         {main_output_dir}/multiqc_report.html
-
+        
         echo "MultiQC complete: $(date)" >> {log}
         """
-    
+
 # Final clean up superfluous files
 rule cleanup_files:
     input:
@@ -3234,9 +3236,6 @@ rule cleanup_files:
             f.write("\nPreprocessing mode: concat")
 
 
-
-
-
 # ----- SE MitoGeneExtractor (mirrors MitoGeneExtractor_concat) -----
 rule MitoGeneExtractor_se:
     input:
@@ -3255,7 +3254,6 @@ rule MitoGeneExtractor_se:
         t=t,
         output_dir=barcode_recovery_dir_se,
         vulgar_dir=lambda wildcards: os.path.join(barcode_recovery_dir_se, f"logs/mge/{wildcards.sample}_r_{wildcards.r}_s_{wildcards.s}/")
-    threads: rule_resources["MitoGeneExtractor"]["threads"]
     resources:
         mem_mb=lambda wildcards, attempt: rule_resources["MitoGeneExtractor"]["mem_mb"] * attempt,
         slurm_partition=rule_resources["MitoGeneExtractor"]["partition"]