feat: add automatic PlusPF re-screening for high-unclassified samples - #319
Open
adRn-s wants to merge 17 commits into
Open
feat: add automatic PlusPF re-screening for high-unclassified samples#319adRn-s wants to merge 17 commits into
adRn-s wants to merge 17 commits into
Conversation
runPlusPF() now runs kraken2 strictly serially (not the configthreads // 5 pooling formula kraken() uses) and adds --memory-mapping, since the PlusPF index is ~75-80GB and loading several copies concurrently risks OOM-killing the demux server. kraken()'s escalation path now bails out (with a log line, no crash) when plusPFdb is missing, empty, or points at a nonexistent path, and when ssdf has no Library_Type column. The escalation loop no longer unpacks krakenfqs()'s result unconditionally, since it can return None. screening.pickThreshold() uses .get()/.getfloat() with fallbacks instead of bare indexing. screening.parseUnclassifiedPct() now also catches ParserError, KeyError, and OSError, not just EmptyDataError.
…-rail gaps Adds an integration test exercising runPlusPF() writing a report and buildContaminationDic() reading it back through the same *.plusPF.krakenreport filename, rather than each side only asserting against its own hardcoded literal. Also adds coverage for kraken()'s new guards: an ssdf with no Library_Type column, and a plusPFdb pointing at a nonexistent (e.g. template placeholder) path.
The screening heading used ^^^^^ (already claimed by a deeper heading elsewhere in the file), which caused misc/communication to render as its children in the TOC instead of siblings; switched to ---- to match every sibling section. Also: reworded executables.rst's PlusPF paragraph, which incorrectly claimed the escalation result appears in the multiQC report -- it is deliberately excluded via fn_ignore_files and only ever shown in the core-team email. Fixed a stray %% (configparser-escaping leftover) rendering literally in the built HTML. Documented that omitting [screening], or leaving plusPFdb unset/pointing at a nonexistent path, turns escalation off without touching any other config, and noted plusPFdb should be sized to fit the server's available RAM for --memory-mapping to be effective.
krakenfqs() raises IndexError (not None) for a sample folder with zero matching fastq files; the escalation loop's guard only covered the too-many-fastqs None-return case. Also relax runPlusPF()'s pooling from strictly serial to half of kraken()'s configthreads // 5 ratio -- escalation is rare enough in practice (~1-2 flowcells/week) that this keeps a wide safety margin against the PlusPF index's size while no longer forcing every run onto a single worker.
Email-only surfacing turned out to be insufficient in practice. Adds a distinct MultiQC custom_content table (id 'plusPF_escalation') built by the new plusPFEscalationTable(), listing each escalated sample's top PlusPF hit and its percentage. The table only appears on a project with at least one escalated sample, and is written/removed by md5_multiqc() the same way the existing parkour/seqreport/index custom_content files are. This stays a separate section rather than folding it into the routine kraken module's output -- the raw .plusPF.krakenreport is still excluded via fn_ignore_files, so the underlying report can't get auto-detected as a phantom extra sample; the new table is the only thing MultiQC ever sees for it.
clmpRunner previously assumed tmp.fq.gz always existed after clumpify ran, and never checked its exit code. A silent clumpify failure (stderr is discarded) left no tmp.fq.gz, so splitFastq ran against nothing and the final os.remove raised an unhandled FileNotFoundError that killed the whole multiprocessing pool instead of the existing failure path in clumper() (mailHome + exit).
plusPFEscalationTable() only reported each escalated samples single top hit (by direct read count) and its cumulative clade percentage. Replace it with plusPFEscalationBargraph(), which stacks the top-N taxa (summed by direct read count across escalated samples) plus Unclassified and an Other bucket, matching the shape of the routine Kraken2 module plot. Values are computed from each krakenreport rows *direct* read count rather than its cumulative clade percentage: direct counts partition every read exactly once, so the stacked bars sum to ~100% per sample. Cumulative clade percentages overlap between parent and child taxa and would double-count if stacked directly. Also place the plusPF_escalation module right after the routine kraken modules own section via report_section_order, so it reads as a follow-up to it rather than an unrelated section elsewhere in the report.
Replace the single-dataset bargraph (one fixed top-5 list mixing all ranks) with a multi-rank tabbed bar plot -- one tab per taxonomic rank (Species/Genus/Family/.../Domain/Unclassified) -- mirroring the layout of MultiQC built-in Kraken modules own Top taxa plot exactly, down to the same Percentages/Counts toggle and Unclassified/Other categories. A plain custom_content TSV can only carry one dataset, so it cannot reproduce those rank-switching tabs. Switch plusPFEscalationBargraph() to emit a custom_content JSON payload instead (list of datasets + list of categories + pconfig.data_labels), which is what the built-in module itself is built on. postmux.py now writes/removes a .json file for this instead of a .tsv. Top-N taxa are now picked per rank (matching the native modules own approach) rather than mixing ranks together in one flat top-5, which was hiding higher-rank signal (e.g. family-level Enterobacteriaceae, which often gets far more direct reads than any single species below it once reads cant be resolved further).
Test_plusPFEscalationTable still targeted the removed plusPFEscalationTable function, breaking CI collection since it was replaced by plusPFEscalationBargraph in 31bb998/329e4ed.
…n-screening # Conflicts: # tests/test_misc.py
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
1 participant
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.
Pull Request Names
Follow Conventional Commits
Briefly:
Make sure the PR has a title that adheres to the following scheme:
If functionality is broken, append a
!to the type, e.g.fix!:orfeat!:, to indicate a breaking change.Description of Changes
Adds automatic escalation screening: when a sample's routine kraken2 screen (against the small curated
contaminomedb) reports more unclassified reads than a configurable threshold, dissectBCL re-screens just that sample against a much broader Kraken2 "PlusPF" index, and surfaces the result in the core-team email report.[screening]config section (dissectBCL.ini, documented indocs/config.rst):plusPFdb,unclassified_threshold(default 10%),relaxed_library_types(defaultATAC-Seq),relaxed_threshold(default 20%). The threshold values were derived empirically from ~13 months of real flowcell QC data (percentile analysis + per-protocol breakdown) — see commit history / PR discussion for the full derivation.src/dissectBCL/screening.py: pure logic to parse a kraken2 report's unclassified %, pick the threshold for a sample'sLibrary_Type, and decide whether to escalate.src/dissectBCL/postmux.py: newrunPlusPF()runs the escalation kraken2 pass (serially, with--memory-mapping, given the PlusPF index's size);kraken()now callsscreening.needsEscalation()per sample and triggersrunPlusPF()for flagged ones. Degrades to a no-op (never aborts a flowcell) when[screening]is absent,plusPFdbisn't configured/doesn't exist, orLibrary_Typeisn't available.src/dissectBCL/fakeNews.py/flowcell.py: escalation results surface as aplusPFcolumn in the core-team email table. Deliberately excluded from the multiQC report (fn_ignore_filesinsrc/dissectBCL/misc.py) so the escalation report never appears as a phantom extra sample.plusPFdbindissectBCL.iniis pointed at a real, downloaded PlusPF kraken2 index (see https://benlangmead.github.io/aws-indexes/k2) — until then it's a safe no-op.Test plan
pytest— full suite passing (89/89)ruff check/ruff format --checkcleansphinx -b html), no new warningsNone-unpack edge case, and a docs/implementation mismatch)plusPFdbstill needs to be deployed there before this can be exercised end-to-end in production