Skip gene_summary header in mageckGSEA (fixes phantom "id" gene) - #21
Merged
Conversation
… gene
mageck2 pathway --method gsea passes a *.gene_summary.txt whose first row is a
header, but mageckGSEA's parserankfile read every line, so the header became a
phantom gene named "id" (score atof("neg|score") == 0). That inflated the gene
count by one and shifted the rank transform, percentile, and permutation math
for every pathway on the default workflow. The RRA path already skips the header
(mageck_readgeneranking), so only the C++ GSEA path was affected.
Add a -H / --skip-header switch to mageckGSEA that drops the first line, and pass
it from all three gseacommand builders in pathwayFunc.py (the gene-ranking input
always carries a header, matching the RRA reader's assumption).
Adds a regression test asserting a headered file read with -H scores identically
to the same data with no header line, plus a CHANGELOG entry. Addresses the P2
Codex review comment on PR #18.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EZcdrQdWJfczuqaMeZ3DBD
The prior assertion checked ES != 0.0 on a top-3-genes pathway whose enrichment score is ~1.4e-17 (floating-point zero). macOS prints it as 1.38778e-17 while Linux prints 0.0, so the sanity check failed on Ubuntu CI even though the real assertion (headered + -H == headerless) passed on every platform. Use a pathway (G14..G18, the lowest-scored genes, which cluster at one end of the ranking) that yields a robustly non-zero ES (~0.87), and assert ES > 0.5. The core regression -- -H makes a headered file score identically to the same data headerless -- is unchanged. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EZcdrQdWJfczuqaMeZ3DBD
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.
Why
Addresses the P2 Codex review comment on #18: "Skip the gene_summary header before ranking."
mageck2 pathway --method gsea(the default) passes a*.gene_summary.txtwhose first row is a header (id,num,neg|score, …).mageckGSEA'sparserankfile()read every line, so the header became a phantom gene namedidwith scoreatof("neg|score") == 0. That inflated the total gene count by one and shifted the rank transform, percentile, and permutation calculations for every pathway on the default workflow — a small but systematic bias.The RRA path is unaffected: its Python reader (
mageck_readgeneranking) already skips line 1. Only the C++ GSEA path had the bug.What
-H/--skip-headerswitch tomageckGSEA;parserankfile()drops the first line when set.-Hfrom all threegseacommandbuilders inpathwayFunc.py. The gene-ranking input always carries a header (consistent with the RRA reader), so it's passed unconditionally.-Hscores identically to the same data with no header line.Verification
Clean-venv source install compiles RRA + mageckGSEA onto PATH; with
-H, a headered 4-gene file reportsrank size:4(header dropped, previously 5). Full smoke suite 14/14 passes.Note on the other Codex comment
The P1 comment on #18 (
operator[]mapping off-list genes to rank 0) was already fixed in #19 and is onmainvia #20 — Codex reviewed the pre-#19 commit. No action needed there.🤖 Generated with Claude Code