Repository navigation
v3.0.7 test branch (new_human_filter) merge - #14
Merged
Merged
Conversation
01_human_cox1_filter.py could not detect human reads. It compared read[:n] against HUMAN_COX1[:n] with both anchored at index 0, with no alignment or offset search, after degapping the MGE-aligned read and so discarding the only positional information available. A verbatim 300bp slice of its own reference scored 1.000 at offset 0 and 0.19-0.31 at every other offset, against a random baseline of 0.249 and a threshold of 0.95. Only a read beginning at exactly COX1 base 1 could ever be removed, so cleaning_removed_human was 0 on every sample and human consensuses passed through to taxonomic validation. Replace it with 01_human_mitogenome_filter.py, which maps reads to a human reference with bwa aln. Measured on a heavily contaminated sample, mapping to a mitogenome + NUMT reference removes 99% of human reads while keeping every target read; the previous mitogenome-only k-mer approaches removed 85%. The reference is a config value (fasta_cleaner.human_reference) rather than a committed file, defaulting to the packaged human mitogenome. A mitogenome-only reference cannot catch NUMT-derived reads by construction, so the config documents pointing it at a mitogenome + NUMT FASTA. fasta_cleaner.human_threshold is removed: mapping has no similarity threshold. Indexing moves out of the script into the bwa_index_human_ref rule, which runs once per run and is shared by the merge, concat and se filter rules. Building on demand inside the script would have re-indexed per invocation and let concurrent mode jobs race to write the same index files. The index is written to resources/bwa_index/ in the working directory via `bwa index -p`, never beside the reference, which is read-only in a pip-installed package; bwa aln, samse and mem read only the index files and never reopen the FASTA, so no copy or symlink of the reference is needed. Outputs are temp() because the default reference re-indexes in about a second; keep_bwa_index: true covers genome-scale references, where bwa index costs an hour or more per gigabase. The metrics CSV keeps the column order and FILE_SUMMARY sentinel that 06_aggregate_filter_metrics.py:parse_human_metrics() requires. That parse reads only rows whose sequence_id is FILE_SUMMARY and keys samples off base_name, and is wrapped in a bare try/except, so a schema change there reports removed_human = 0 for every sample with no warning - the same symptom this change exists to fix. METRICS_COLUMNS and build_summary_row() make the contract explicit and testable. Output paths, file names and metric names are unchanged, so the all targets, at_content_filter rules and metrics aggregation are unaffected. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013enV9Cqh1egXMhrpRtDyCP
nhmmer reports one envelope per contiguous stretch of the query matching the profile, so a consensus with a low-coverage hole in the middle comes back as two envelopes. Both are real data. parse_nhmmer_result() held a single best_result and overwrote it only on a strictly lower E-value, so every other envelope was logged as "Found alignment" and then dropped with no further log line at all - only insignificant envelopes got an explicit "Rejected" line, so a silently discarded significant envelope looked accepted in the log. The effect was to truncate the barcode to whichever fragment happened to score best. On UK016-H10 r_1.5_s_50_fcleaner_concat, a 299-column envelope at E=6.1e-63 lost to a 291-column one at E=3.9e-73 and the pipeline emitted 290 bases from a consensus whose longest unambiguous stretch was 900. Eight rows for that sample reported exactly 290 bases - all the same surviving envelope, HMM 367-657. Barcode selection was therefore driven by how nhmmer happened to fragment the match rather than by sequence quality: the row that won selection did so because its match came back as one envelope instead of two. parse_nhmmer_result() now returns every significant envelope, chosen greedily in E-value order and accepted only where it overlaps nothing already accepted in BOTH HMM and sequence coordinates. Both tests are needed: overlapping HMM coordinates would let one envelope silently overwrite another's model columns, and overlapping sequence coordinates would place the same query bases at two different model positions. Coordinates are normalised before comparison so a minus-strand hit, which nhmmer reports with seq_from > seq_to, is not missed. Because the best envelope is always accepted first, the result can never cover less of the model than before - any barcode that gets shorter is a bug, not a trade-off. construct_hmm_space_from_alignment() places each envelope at its own model coordinates, so the span between them stays a gap, becomes internal N via replace_gaps_with_n(), and counts against barcode_base_count rather than being mistaken for called sequence. The reading frame survives that span because HMM positions are profile columns. Envelopes dropped for overlap are now logged, naming the envelope they clashed with and the coordinate system involved. None rather than [] is still returned when nothing is significant, preserving the contract run_nhmmer_on_sequence()'s callers test with `is None`. Verified against the real UK016-H10 tabular rows: r_1.5_s_50 recovers HMM 2-657 (656 positions, 66 N between envelopes) in place of 290; r_1.5_s_100 recovers HMM 119-657 in place of 290; r_1.3_s_100 and the uncleaned r_1_s_50_concat, which each matched as a single envelope, are unchanged. This changes barcode_length, barcode_base_count and barcode_rank for any sample whose match was fragmented, and can change which barcode is selected. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013enV9Cqh1egXMhrpRtDyCP
The envelope work in 8c44f22 merges every significant, non-overlapping envelope. On the UK016 plate that rescued E11 and G07 but lost C08, which had passed before. C08 r_1.3_s_50_fcleaner_concat shows the failure mode. nhmmer reported two significant envelopes: HMM 233-619 (387 positions, E=2e-83, score 269.1) and HMM 2-50 (49 positions, E=3.6e-06, score 14.3). They are 182 model positions apart, and that gap has no data, so merging filled it with N. The barcode grew 231 columns but only 48 of them were informative: 387 + 182 + 49 = 618 length 18 + 182 + 1 = 201 N 369 + 48 = 417 bases passes_quality_criteria() rejects a barcode whose ambiguous/called ratio reaches 0.30, and the ratio went from 0.049 to 0.482. A sample that had produced a 99.251% match to Triphleba minuta was discarded in favour of a longer, emptier sequence. An extra envelope is now merged only when its model span exceeds the gap it would open. An envelope lying inside the accepted span overlaps nothing (those are rejected earlier) and fills a hole that is already N, so it costs nothing and is always taken. The best-E-value envelope remains the anchor, so no barcode can come out shorter than before the envelope work; a shorter barcode would be a bug, not a trade-off. Score and bias were evaluated as alternatives and rejected. Both describe an envelope in isolation and are blind to where it sits relative to the rest of the barcode, which is what sets the N cost: a 49-position envelope 5 positions away and one 182 positions away have identical score and bias and opposite consequences. Measured over the 430 secondary envelopes in the plate log, the best single score threshold reproduces the net-gain verdict on 77% of them and the best bias/score threshold on 68%; the two classes overlap completely (dropped envelopes reach score 103.6, kept ones fall to 12.4). Raw bias selects backwards, since it scales with length - G07's kept HMM 3-318 has bias 39.1 while C08's spurious HMM 2-50 has 3.4. Verified by replaying all 2,242 sequence blocks from the plate's structural_validation.log through the new parser: 214 rows change envelope set, no sample loses every row that passes the quality gate, and the documented cases come out exactly as intended - C08 reverts to [(233,619)] while G07 [(3,318),(373,657)] and E11 [(2,263),(361,627)] are untouched. Both branches log the span-versus-gap comparison so the decision is auditable. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013enV9Cqh1egXMhrpRtDyCP
…refactor Snakefile, and streamline comments in configs
Updated CHANGELOG to reflect changes in version 3.0.7
…-inclusive default in README Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011p7WP8g96BFFSiyAtAWBAZ
This document provides a comprehensive reference for human contamination, detailing the construction of a human reference FASTA file that includes both the human mitogenome and NUMTs. It also outlines the process for using a different reference in the configuration.
…ngelog publish.yml runs on a published GitHub Release: reuses ci.yml for tests, checks the tag matches the pyproject version, builds and twine-checks the sdist and wheel, smoke-tests the installed wheel, then uploads to PyPI via Trusted Publishing. The smoke test exposed that the CLI imports yaml while the package declared no dependencies, so pyyaml is now declared. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011p7WP8g96BFFSiyAtAWBAZ
SchistoDan
changed the title
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.
structural_validation.pynow 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.01_human_cox1_filter.pywith01_human_mitogenome_filter.py, to remove human contamination by mapping reads to a human mitogenome reference and NUMTs withbwa alninstead of comparing them positionally against a hardcoded COX1 sequence.fasta_cleaner.human_referenceto set the mapping reference.nulluses the packaged human reference.fasta_cleaner.human_thresholdis removed - mapping has no similarity threshold.bwa_index_human_refrule, which indexes the human reference once per run and is shared by the merge, concat and se filter rules. The index is written toresources/bwa_index/in the working directory viabwa index -p. It is removed at the end of the run; setkeep_bwa_index: truewhen using a genome-scale/alternative reference.human_cox1_filterrules and therules.human_cox1_filterresource key tohuman_mitogenome_filter. Output paths, file names and metric names are unchanged.show_analysis_paths: falseinmultiqc_config.yaml), so server directory layouts (e.g. Galaxy job directories) are not exposed to end users. Issue #13filecolumn ofstructural_validation.csvnow holds paths relative to the output directory (e.g.03_barcode_recovery/concat_mode/consensus/...) instead of absolute paths. Issue #13val_barcode_recovery_metrics_mergerule toval_barcode_recovery_metrics_combineto be more descriptive (not only a product of the 'merge' pre-processing mode).gene_fetch.typeis now checked at startup with the other requiredgene_fetchkeys, so a missing value fails immediately with a clear message instead of a KeyError when the rule is built.pyyamlas a package dependency..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.