Merge BLAST chunking, fix temp-file leak and summary-CSV bug, and improved rule resources - #12
Merged
Merged
Conversation
tv_local_blast.py ran one blastn process per query sequence, so a large reference database was loaded and discarded once per barcode. For a 500-barcode run against BOLDistilled that fixed startup cost dominates the search itself. Queries are now batched into chunks of CHUNK_SIZE (50) sequences per blastn invocation, with up to --processes chunks running concurrently; the final chunk holds the remainder. A 120-sequence run drops from 120 invocations to 3. Chunks are written with synthetic query IDs (>q1, >q2, ...) so results demultiplex back to their source sequence by exact lookup rather than by re-deriving the name from BLAST's qseqid, which is only the first whitespace token of the header. The real query ID is restored in column 0 before each per-sequence TSV is written, so output is unchanged. Verified byte-for-byte identical per-sequence TSVs and summary CSV against the previous implementation over 120 sequences, including sort stability across equal-percent-identity ties (which downstream taxonomy selection depends on). Other changes in the same path: - Resume moves out of the worker into the input scan, so finished sequences no longer cost a temp file and a pool task. Per-sequence TSVs are written via a .part file and renamed, so an interrupted run cannot leave a truncated file that a later run treats as complete. - ThreadPoolExecutor replaces ProcessPoolExecutor. The work is entirely subprocess-bound, and this removes the per-task pickling of the growing processed_sequences dict. That dict is now populated in the main thread; the previous in-worker mutation was silently discarded in the child process. - A failed chunk retries its sequences individually so one bad record cannot lose its 49 neighbours. Sequences that still fail get no output file and cause a non-zero exit with no summary CSV, rather than being handed to tv_blast2taxonomy as a genuine no-match. Successful TSVs stay on disk, so retries converge. --allow-partial-failures restores the previous leniency. - Fixed the summary CSV dropping every hit for sequences whose header carries a description: it is keyed on the sanitized full header but was looked up by qseqid, so 'SEQ1 some description' never matched 'SEQ1'. - Removed the unused blast_options entry from the Snakefile's taxonomic_validation fallback and noted that the fallback is unreachable. - Corrected the documented BLAST hit limit (100, not 500) and the config.yaml threads comment, which described 4 threads per BLASTn process where the script uses -num_threads 1. Adds unit tests for chunk boundaries, the input scan, demux, tie stability and the failure path, plus a blast_runner fixture that bypasses the blastn version check so the suite still needs no BLAST binary. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013j5MNgaVXKBo7QaqBVhVqZ
…ight-size resources
run_nhmmer_on_sequence created its per-sequence query FASTA and nhmmer tabular
output with tempfile.NamedTemporaryFile(delete=False), and the script never
unlinked either. Two files per sequence were left in TMPDIR for the lifetime of
the job. That is 48 per sample, since every sample is validated across 6
MitoGeneExtractor parameter combinations and 4 input FASTAs (pre- and
post-clean consensus for both concat and merge mode) - roughly 24,000 files for
a 500-sample run. SLURM usually wipes a per-job TMPDIR, which is why this has
gone unnoticed, but it can fill node-local /tmp on a large run.
Both files now live in a single tempfile.TemporaryDirectory() per call, so they
are removed on exit including on error. An absent tabular file is treated as
"no hits" rather than falling through to the FileNotFoundError handler, which
would have misreported it as a missing nhmmer binary.
Adds --threads (default 1), plumbed through align_sequence_with_nhmmer and
analyse_fasta to nhmmer's --cpu, and passed as --threads {threads} from the
rule. The default keeps standalone behaviour identical. Note this sets threads
per nhmmer invocation, not concurrent invocations - the script still processes
sequences serially, one nhmmer call each.
Right-sizes the rule's resources from 32 GB / 1 thread / himem to 4 GB /
8 threads / medium. Peak memory scales with sequence count, not read depth:
every result retains two Biopython SeqRecords until clean_results_for_csv_output
runs at the very end, measured at 4.67 KB per sequence, so roughly 110 MB at
1,000 samples and 550 MB at 5,000. The himem partition and the 32 GB request
were over-provisioned by about two orders of magnitude.
Adds five tests covering the --cpu wiring and temp-file cleanup on both the
success and failure paths. Three of them fail against the previous
implementation, which is what makes them regression tests; they use a fake
subprocess.run so CI still needs no nhmmer binary.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013j5MNgaVXKBo7QaqBVhVqZ
… one
Audited every entry in the rules: block against its Snakefile rule and the script
or tool it invokes, then rewrote the block with per-rule guidance.
Two structural facts now recorded in the file:
Only fastp_qc, clean_headers_merge, fastq_concat, quality_trim, downsample and
MitoGeneExtractor scatter over samples. Every other rule - including all the
fasta_cleaner filters - is a single job for the whole run. The block is split
into PER-SAMPLE and AGGREGATE sections accordingly, because the first group costs
the same for a 1-sample run as for a 500-sample one and only the second needs
headroom as runs grow.
Eleven rules reserved CPUs their commands never consume: gene_fetch (NCBI
API-bound), downsample (reformat.sh is invoked without a thread parameter),
gzip_merged_clean (a serial gzip loop in a run: block), blast2taxonomy,
extract_stats_to_csv, download_taxdump, multiqc_plots, multiqc, clean_headers_merge
and fastq_concat (no {threads} in the shell), and MitoGeneExtractor (Exonerate is
single-threaded; {threads} appears only in an echo). Those drop to 1. Conversely
reference_filter was left at 1 despite having the same process pool as the other
filters, so it ran serially for no reason; it goes to 4.
Memory figures are measured rather than assumed where they scale with run size:
blast2taxonomy holds every hit for every barcode at ~44 KB per barcode (~340 MB at
500 samples), and structural_validation at ~110 MB per 1000 samples. The
fasta_cleaner filters parallelise over per-sample files, so their peak is
threads x largest single alignment - it tracks read depth, not sample count, which
is the guideline most likely to be misread.
Net effect is 249 GB -> 145 GB of reserved memory and 73 -> 55 threads, with
nothing raised except reference_filter. The header explains that mem_mb * attempt
retry scaling is the intended flexibility mechanism, so base values are sized for a
typical run rather than the worst imaginable case, and warns that the four
partition keys the Snakefile indexes directly must not be deleted.
Not changed, but flagged in passing: downsample sets no -Xmx, so the BBTools JVM
sizes its heap from visible node RAM rather than the cgroup and can overshoot its
request; and MitoGeneExtractor inherits runtime=720 from the SLURM profile, which
is a profile setting rather than a rules: one.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013j5MNgaVXKBo7QaqBVhVqZ
config.yaml deliberately no longer declares threads for the rules whose commands never consume it, but test_rules_have_threads asserted every rule declares one, so it enforced the opposite of the intent and CI was red. It is replaced by test_rules_threads_valid_when_present, which checks the value is a positive int only where the key exists. test_rules_have_mem_mb is unchanged - mem_mb really is required by every rule. Also adds test_snakefile_resource_reads_resolve, covering the failure mode that removing those keys actually caused: MitoGeneExtractor_se and downsample_se kept reading rule_resources[...]["threads"] after the config entries were deleted. Snakemake evaluates those reads at parse time, so the MitoGeneExtractor_se one (defined at column 0, unconditional) raised a KeyError on every run, PE or SE, before any job was scheduled. Nothing in the suite caught it: the Snakemake dryrun test self-skips when snakemake is absent, which is the normal state in CI. The new test parses the Snakefile for rule_resources["<rule>"]["<key>"] and asserts each resolves against config.yaml. It needs no snakemake and reports the offending line numbers. Verified against the pre-fix Snakefile, where it fails naming both Snakefile:1188 and Snakefile:3257. The rule-name character class is [A-Za-z_] rather than [a-z_] on purpose: MitoGeneExtractor is the only rule with uppercase in its name, and a lowercase-only pattern silently skips exactly the read that broke every run. That is asserted in the docstring so the next person does not "tidy" it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013j5MNgaVXKBo7QaqBVhVqZ
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary