Resolve dropped duplicate IDs when matching --pg-pair-only - #31
Conversation
Counting resolves a guide sequence to the one surviving ID for that sequence, but --pg-pair-only was matched against the ID exactly as written in the pair file. A pair naming an ID that had been dropped for duplicating an earlier sequence could therefore never match, and the whole construct was excluded from pg_count.txt -- even though its sequence pair had been counted correctly all along. For a dual-guide library indexed by construct rather than by distinct guide, that is most of the file. The remedy was to rewrite the pair file by hand against the surviving IDs, which issue #15 shows is mechanical. Do the remapping instead, using the {dropped_id: representative_id} map the loader now produces. It is carried on args so the counting stage matches the same names it labels rows with. Report it: the per-line warning now says the ID was resolved rather than that the pair will be excluded, and a summary states how many pair-file entries were remapped and what they became. Resolution can collide. If both libraries contain duplicated sequences, two distinct pair-file entries can resolve onto one row, and their reads really are indistinguishable -- the sequence pair is the same. That is a property of the data, not something to fix, but it must not happen silently, so name both entries and say their reads are reported as a single row. Tested end to end against a construct-indexed library: the pair naming a dropped ID is now counted under the representative ID, the summary is logged, and the colliding case warns while still producing one merged row. Closes #28 Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
https://github.com/davidliwei/mageck2/blob/3d852eea699432ea0e28ec4cc2e3478bc48748ea/mageckCount.py#L384-L386
When a construct-indexed pair file contains many unique dropped IDs—the large-library case this change is intended to support—each alias performs a linear membership scan over every previously collected alias and then extends the list, making pair-file loading O(R²). With tens or hundreds of thousands of remapped records, writing pg_count.txt can become prohibitively slow; track uniqueness in a set while retaining only the five examples needed for logging.
https://github.com/davidliwei/mageck2/blob/3d852eea699432ea0e28ec4cc2e3478bc48748ea/mageckCount.py#L381
If an ID is first dropped because its sequence duplicates an earlier row and a later row reuses that ID with a distinct sequence, mageckcount_checklists accepts the later row but leaves the ID in the dropped map. This unconditional lookup then redirects a valid surviving ID to the earlier representative, causing the pair file to allow the wrong sequence pair and filter the intended one; aliases should only be applied when the ID is absent from the final loaded library, or shadowed entries should be removed from the dropped map.
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
Two problems in the resolution added on this branch, both raised in review. An ID dropped for duplicating an earlier sequence can be admitted by a later row carrying a sequence of its own: the loader's duplicate-ID guard tests what has already loaded, and the dropped ID is not in it. That ID is then a real sgRNA while still holding a stale entry in the dropped map, so resolving pair-file IDs through the map unconditionally redirected a valid ID to an unrelated representative -- allowing the wrong sequence pair and filtering the intended one, silently and with plausible-looking counts. Remove such entries from the map once loading is done: an alias is only meaningful for an ID the library does not have. Collecting the aliases for the log line scanned a growing list for each record, which is quadratic in the number of remapped entries -- worst on exactly the construct-indexed libraries this feature exists to support. Track uniqueness in a set and retain only the handful of examples the message prints. Reporting is unchanged; at 20k unique aliases the inner loop goes from 1.17 s to 2.5 ms. Refs #28 Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
Both findings are correct and both are fixed in 5037752. P2 — remapping IDs that survive loading. Confirmed. Fixed at the source rather than the use site: entries whose ID is present in the final library are removed from the map once loading finishes, so an alias only ever exists for an ID the library does not have. That also keeps the duplicate-drop report honest — it no longer names an sgRNA that survives. Regression test P1 — quadratic alias collection. Confirmed: Now a set for uniqueness plus at most five retained example strings, as suggested. Reporting semantics are unchanged — the trailing Suite is 31 passed, up from 30. |
|
@codex review |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 503775270e
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| continue | ||
| nrecord+=1 | ||
| sg_1=field[0]; sg_2=field[1] | ||
| sg_1=dropped.get(field[0],field[0]); sg_2=dropped2.get(field[1],field[1]) |
There was a problem hiding this comment.
When a library repeats the same dropped sgRNA ID for two different sequences that are each represented by earlier rows, mageckcount_checklists overwrites dupmap[id] on every duplicate-sequence row. This lookup therefore resolves the pair-file ID to whichever representative appeared last, causing that representative's sequence pairs to pass --pg-pair-only arbitrarily and making results depend on library row order. Track multiple representatives for a dropped ID and reject or leave ambiguous IDs unresolved instead of selecting the last one.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Correct, and fixed in 7cb6496.
Reproduced before changing anything: with sgA/sgB claiming two sequences and sgX duplicating both, the map came out as {'sgX': 'sgB'} — and reversing sgX's two rows gives sgA. Neither row loads sgX into the library, so the previous commit's survivor cleanup does not catch it, and each duplicate-sequence row overwrote the entry.
There is no right answer for such an ID: it names two different sequences and nothing in the data prefers one. So rather than picking deterministically, the competing representatives are recorded, the ID is left out of the map, and a warning names each ID with the representatives that compete for it:
1 sgRNA IDs in lib.txt each name more than one duplicated sequence, so no
single sgRNA represents them. They are left unresolved, and pairs naming
them will not match --pg-pair-only: sgX -> sgA/sgB
A pair file naming such an ID will not match, which is the honest outcome, and it is now stated at load time instead of surfacing later as a mislabelled count. Ambiguous IDs are removed before the duplicate-drop summary is emitted, so that message no longer claims a representative it does not have.
test_ambiguous_duplicate_id_is_left_unresolved builds the library in both row orders and asserts the same result for each, since order dependence was the defect.
Suite is 32 passed, up from 31.
An sgRNA ID can appear on more than one duplicate-sequence row, each naming a different sequence that an earlier row already represents. Every such row overwrote the ID's entry in the dropped map, so it resolved to whichever representative the library happened to list last -- and which sequence pair then passed --pg-pair-only depended on library row order. There is no right answer for such an ID: it names two different sequences, and nothing in the data prefers one. Record the competing representatives, leave the ID out of the map, and warn, naming each ID and the representatives that compete for it. A pair file naming it will not match, which is the honest outcome and is now stated up front rather than discovered as a silently mislabelled count. The test asserts the same result for both row orders, since order dependence was the defect. Refs #28 Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
@codex review |
|
Codex Review: Didn't find any major issues. Swish! Reviewed commit: ℹ️ About Codex in GitHubYour team has set up Codex to review pull requests in this repo. Reviews are triggered when you
If Codex has suggestions, it will comment; otherwise it will react with 👍. Codex can also answer questions or update the PR. Try commenting "@codex address that feedback". |
1 participant
Closes #28. Builds on the
{dropped_id: representative_id}map added in #30.The problem
Counting resolves a guide sequence to the one surviving ID for that sequence, but
--pg-pair-onlywas matched against the ID exactly as written in the pair file. A pair naming an ID dropped for duplicating an earlier sequence could never match, so the whole construct was excluded frompg_count.txt— even though its sequence pair had been counted correctly all along. On a dual-guide library indexed by construct rather than by distinct guide, that is most of the file.What changed
Resolution. Pair-file IDs are resolved through the dropped-ID map before building the allowed-pair set, so
--pg-pair-onlymatches the same namespg_count.txtlabels rows with. The map is carried onargsfrommageckcount_checkargsto the counting stage.Reporting. The per-line warning from #30 said the pair would be excluded and to edit the pair file; both are now obsolete, so it says the ID was resolved and the pair is counted under it. A summary states the scale:
Collisions. If both libraries contain duplicated sequences, two distinct pair-file entries can resolve onto one row. Their reads genuinely are indistinguishable — the sequence pair is identical — so merging them is correct, but silence is not:
Help text.
--pg-pair-onlynow states that dropped IDs are resolved, so a pair file written against the original nomenclature does not need rewriting.Tests
Four new tests, each watched fail before the code that passes it, on a fixture shaped like the library in #15 (
sg1andsg1_dupsharing a sequence, pair file namingsg1_dup):Full suite: 30 passed (26 on
main).Note on scope
This does not restore the original construct IDs as row labels — the
sequence -> [ids]multimap discussed in #15 and recorded under "Not doing" in #28. Counts are already correct; only the label differs, and that is a downstream rename. The reasoning is unchanged.Follow-up
#29 still stands: extend the paired-guide fixture to
--reverse-complement-2,--pairguide auto, and--pg-min-read, and fix the unboundedctab_umi[None]accumulation insecondpairmode.