v2.3.0 PR - #116
Conversation
…y slow to account for an extremely rare failcase
There was a problem hiding this comment.
The FLUSH_BYTES bit of this makes the whole thing read a little more convoluted, but batching writes makes sense
There was a problem hiding this comment.
Yeah it's fundamentally a better way to handle the problem imo but the implementation could probably do with tidying up...
|
Changing rank -> simple rank: do the sankeys in the html report make sense still (I haven't run the code while reading it, and I'm not sure anyone looks at them but want to know anyway...). In our first iteration of scylla we were extracting all things that were at simple rank S (including S, S1, S2 etc) as a way of capturing viruses and bacteria where the named thing wasn't necessarily the species, but I remember there being some strangeness in the sankeys, and also there being SO MANY almost identical extracted read files. |
|
Just out of interest, how much does rammap improve the performance of HCID detection module? Why is the rust implementation better? |
rmcolq
left a comment
There was a problem hiding this comment.
Caveat, I've not run the code so you may want a second review. This looks good to me. Lots of things combined into one PR, but a lot of them were touched by the replacement of 'copy' output file mode. The code to do with extracting reads is still probably the most fragile bit of the code base and extensively changed, but as all the tests pass and you've checked the outputs match exactly, happy with all the changes.
I don't have exact numbers and this is only really relevant for super large (my testdata was 82G of zymo) but it was fairly significant, 20% faster as a gut estimation. It isn't better in any way other than speed, it produces byte identical outputs. |
Want to double check the extract taxa on some more examples before approving
rmcolq
left a comment
There was a problem hiding this comment.
I'm now happy with the extract code
Reduces runtime and intermediate storage for large inputs, and fixes several taxonomy-parsing bugs surfaced by recent taxonomy versions.
Single-pass read extraction
Previously each output group meant its own process and its own full decompress and parse of the FASTQ. Now:
extract_taxa: one process per sample covering all kreport splits, driven by a--report_configJSON (was one process per split).extract_fractions: 6 processes collapsed to 1.virus_and_unclassified,virusandhuman_filteredare extracted in one pass via--fraction_config.KrakenAssignments.get_read_maps()builds N independent read maps in a single pass over the assignment file.extract_utils: newTaxonWriter(LRU-capped filehandles plus buffered writes) bounds memory and fd use; per-taxon running sums replace per-read quality/length lists; fastermean_phred.HCID
check_hcid.pynow takes read sequences from the SAM rather than re-decompressing the FASTQ, reverse-complementing flag-16 records to restore original orientation. The FASTQ is no longer staged into the process at all.rammap_hcid), with a newdockerfiles/rammap/Dockerfilepulling the pinned v1.1.2 release binary. Verified byte-identical SAM output vs minimap2 2.26.Compression and I/O
bgzip/pigzreplaced withcrabzthroughout (faster, lower CPU, better ratio).fastpnow streams via FIFO/stdout straight into crabz instead of writing an uncompressed temp FASTQ.get_total_length.pyand its two processes, an entire extra FASTQ pass. Total length now comes from the fastcat per-read stats and from the extraction pass itself.check_reads.py: bounded-window duplicate/interleave detection with early exit, plus a long-read fast path.concatenate_reads.py: batched writes; fixedgzip.open(..., mode="w")to"wt".Taxonomy and report correctness
report.py: domain detection by name (modern taxonomy usesR1/R2rather thanD); newsimple_rankso numbered ranks (G1,D2, etc.) match correctly, previously these were missed in per-taxa read files.split_kraken_report.py: ancestor backfill ordering fix.taxonomy.py:childrenas list rather than set;load_entries_from_nodesnow honours itstaxon_idsfilter instead of loading all ~2.6M nodes."Select domain" dropdown showed random mix of taxa for modern taxonomy
NCBI retired
superkingdomin 2024: Bacteria/Archaea/Eukaryota are now rankdomain, Viruses isacellular root, withrealmclades beneath it.aggregate_lineages_bracken.py's hardcodedRANKSallowlist silently dropped any node whose rank wasn't listed, so the real domain nodes vanished and the next allowed rank —kingdom— became the top level of the lineage tree. The dropdown is populated from that tree's top-level keys, so it offered kingdoms and clades: on the test fixture,Pseudomonadati,Opisthokonta(a clade),Heunggongvirae. NowBacteria,Eukaryota,Viruses,Unclassified.RANKSgainsacellular root,domainandrealm, and keepssuperkingdomso pre-2024 dumps still resolve.cellular rootis deliberately excluded — including it would collapse all three cellular domains under a single meaningless node.domain, so it previously contributed to no top-level node at all. Sankey percentages will shift slightly, and the top-level keys in*.kraken.jsonare renamed.bin/ranks.pyholds the vocabulary once (RANKS,DOMAIN_RANKS,UNCLASSIFIED_RANK,KNOWN_DOMAIN_NAMES), imported byaggregate_lineages_bracken.pyandreport.pyand injected into the report's JS bymake_report.py. Three copies previously had to agree by hand, and the JS copy drives both the#rank-selectoptions and the sankey depth filter — where an unrecognised rank yieldsindexOf() == -1and passes every filter, so drift showed up as nodes ignoring the rank selector.sankey.js: nodedomain_namewas derived fromval.rank == "superkingdom", which after the rename matched only the syntheticUnclassifiednode — leaving every tooltip'sDomain:fieldnull. Now tested againstdomainRanks.aggregate_lineages_bracken.py: theUnclassifiednode's rank wasRANKS[0], incidentallysuperkingdom; now an explicitUNCLASSIFIED_RANK.aggregate_lineages_bracken.py:previous_rankwas only assigned on the "node already exists" branch, so thespecies→subspeciesfallback for off-allowlist ranks only fired when the parent species node happened to already be in the tree.report.py:get_domains()re-derived domains fromsimple_rank == "D", which the module's own comment explains is unreliable on modern dumps; it now returns theself.domainsmap already built at load time. Also dropped the strayentry.print()side effect from it and fromget_tips()/get_rank_entries().Containers and config
fastcat/fastp/jqplus crabz;get_versionsreadsconda-metavia jq (Wave images ship no package manager).publish_dir_modeparam (defaultcopy) replacing hardcodedmode: 'copy'.docker-build-push.yml; version bumped from v2.2.1 to v2.3.0.Testing