offtarget: optional per-read XC tags on a target-window BAM for IGV r… - #17
Merged
Merged
Conversation
…eview The ECS caller decides per read whether it supports an edit, aggregates to per-site counts and throws the per-read decision away, so a reviewer sees indel_fraction = 0.42 and has to re-derive by eye which reads made it up. Add --tagged-bam-out to bin/find_edited_reads.py: every read over a target window gets a string tag (default XC) naming its classification, so IGV can colour the pileup by the caller's own reasoning. Two constraints shaped the implementation: * The per-target read loop fetches with a +/-150 bp pad, so overlapping targets make it visit the same alignment record repeatedly. Writing in-loop would emit duplicate records, which IGV renders as inflated depth. Tags are therefore only accumulated in the loop, keyed on (query_name, flag, reference_start) -- query_name alone collides between mates and between primary/supplementary records -- and the BAM is written in a second pass, sorted and indexed. Conflicting calls for one record resolve by a documented precedence (edit > reference > unevaluable > padding). * Output is window-restricted and off by default. On one AAVS1 site14 ECS sample (1149 targets, ~11,000x) it is 800 MB / 21.8 M reads and lifts peak RSS from 1.4 GB to 6.3 GB, so ECS_INDELS gets 24 GB instead of 8 GB while the flag is set. The dominant classification is left implicit rather than stored, which keeps the tag map roughly 40% smaller. Verified on real data: with tagging off the TSV is byte-identical to the unmodified script; with it on, the full-panel BAM has 21,792,054 records, zero missing tags, zero duplicate keys, and is coordinate-sorted and indexable. Note the tag counts do not equal the TSV columns and are not meant to: tags are per alignment record, indel_reads is per fragment and post-filter. Documented. Wires params.offtarget_tagged_bam through ECS_INDELS and GET_INDELS, adds the schema entry, docs with the IGV recipe, and six tests on a new synthetic aligned-read fixture covering dedup, sort order, skip reasons and precedence.
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
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.
…eview
The ECS caller decides per read whether it supports an edit, aggregates to per-site counts and throws the per-read decision away, so a reviewer sees indel_fraction = 0.42 and has to re-derive by eye which reads made it up.
Add --tagged-bam-out to bin/find_edited_reads.py: every read over a target window gets a string tag (default XC) naming its classification, so IGV can colour the pileup by the caller's own reasoning.
Two constraints shaped the implementation:
The per-target read loop fetches with a +/-150 bp pad, so overlapping targets make it visit the same alignment record repeatedly. Writing in-loop would emit duplicate records, which IGV renders as inflated depth. Tags are therefore only accumulated in the loop, keyed on (query_name, flag, reference_start) -- query_name alone collides between mates and between primary/supplementary records -- and the BAM is written in a second pass, sorted and indexed. Conflicting calls for one record resolve by a documented precedence (edit > reference > unevaluable > padding).
Output is window-restricted and off by default. On one AAVS1 site14 ECS sample (1149 targets, ~11,000x) it is 800 MB / 21.8 M reads and lifts peak RSS from 1.4 GB to 6.3 GB, so ECS_INDELS gets 24 GB instead of 8 GB while the flag is set. The dominant classification is left implicit rather than stored, which keeps the tag map roughly 40% smaller.
Verified on real data: with tagging off the TSV is byte-identical to the unmodified script; with it on, the full-panel BAM has 21,792,054 records, zero missing tags, zero duplicate keys, and is coordinate-sorted and indexable.
Note the tag counts do not equal the TSV columns and are not meant to: tags are per alignment record, indel_reads is per fragment and post-filter. Documented.
Wires params.offtarget_tagged_bam through ECS_INDELS and GET_INDELS, adds the schema entry, docs with the IGV recipe, and six tests on a new synthetic aligned-read fixture covering dedup, sort order, skip reasons and precedence.
PR checklist
nf-core lint).nextflow run . -profile test,docker --outdir <OUTDIR>).nextflow run . -profile debug,test,docker --outdir <OUTDIR>).docs/usage.mdis updated.docs/output.mdis updated.CHANGELOG.mdis updated.README.mdis updated (including new tool citations and authors/contributors).