diff --git a/CHANGELOG.md b/CHANGELOG.md index d5c4386..0013a2a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -35,6 +35,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 Previously the permutation test resampled the same full set, found no exceedance of the zero score, and returned `p_permutation=0`; such pathways now short-circuit to a p-value of 1.0 so they cannot distort sorting or FDR. +- `pathway --method gsea` no longer ranks the gene-ranking file's header row as a + phantom gene. The default input is a `*.gene_summary.txt` whose first row is a + header; `mageckGSEA` used to parse it as a gene named `id` with score 0, + inflating the gene count and shifting the ranking, percentile, and permutation + math for every pathway. `mageckGSEA` now accepts a `-H` / `--skip-header` + switch, and `pathway --method gsea` passes it so the header is dropped. ## [0.2.0] - 2026-07-10 diff --git a/gsea/src/calculategsea.cpp b/gsea/src/calculategsea.cpp index 9c0b33b..0a1054c 100644 --- a/gsea/src/calculategsea.cpp +++ b/gsea/src/calculategsea.cpp @@ -21,7 +21,8 @@ struct CallingArgs{ int perms; bool needsort; int scorecolumn; //column for scoring - bool reverse; // whether we need to reverse the order of the gene + bool reverse; // whether we need to reverse the order of the gene + bool skipheader; // skip the first (header) line of the rank file }; typedef map > pathmap; @@ -88,9 +89,11 @@ int parseArguments(int argc, char* argv[],CallingArgs& args){ cmd.add(hassarg); SwitchArg reversearg("e","reverse_value","Reverse the order of the gene."); cmd.add(reversearg); - + SwitchArg skipheaderarg("H","skip_header","Skip the first (header) line of the rank file."); + cmd.add(skipheaderarg); + cmd.parse(argc,argv); - + args.gmtfile=gmtfilearg.getValue(); args.rankfile=rankfilearg.getValue(); args.outputfile=outputfilearg.getValue(); @@ -99,6 +102,7 @@ int parseArguments(int argc, char* argv[],CallingArgs& args){ args.needsort=hassarg.getValue(); args.scorecolumn=scorecolumnarg.getValue(); args.reverse=reversearg.getValue(); + args.skipheader=skipheaderarg.getValue(); }catch(ArgException &e){ cerr<<"error: "<& gname, vector& gscore, Ca return -1; } string oneline; + // Drop the header line if requested (e.g. MAGeCK's gene_summary output, whose + // first row is column names). Otherwise it would be parsed as a phantom gene + // and skew the ranking, percentile, and permutation calculations. + if(args.skipheader){ + getline(ifs,oneline); + } while(true){ getline(ifs,oneline); ncount++; diff --git a/mageck2/pathwayFunc.py b/mageck2/pathwayFunc.py index ad05ce1..91894dd 100644 --- a/mageck2/pathwayFunc.py +++ b/mageck2/pathwayFunc.py @@ -290,6 +290,7 @@ def mageck_pathwaygsa_fast(args): outputfile=args.output_prefix+'.pathway_summary.txt' gseacommand="mageckGSEA " gseacommand+=" -s " + gseacommand+=" -H " # gene_ranking is a gene_summary file with a header row gseacommand+=" -c "+str(columnid)+" " gseacommand+=" -p "+str(args.permutation) gseacommand+=" -g \""+args.gmt_file+"\" " @@ -305,6 +306,7 @@ def mageck_pathwaygsa_fast(args): columnid=args.ranking_column gseacommand="mageckGSEA " gseacommand+=" -s " + gseacommand+=" -H " # gene_ranking is a gene_summary file with a header row gseacommand+=" -c "+str(columnid)+" " gseacommand+=" -p "+str(args.permutation) gseacommand+=" -g \""+args.gmt_file+"\" " @@ -316,6 +318,7 @@ def mageck_pathwaygsa_fast(args): columnid=args.ranking_column_2; # columnid=6 if sgRNA number in positive selection is not omitted gseacommand="mageckGSEA " gseacommand+=" -s " + gseacommand+=" -H " # gene_ranking is a gene_summary file with a header row gseacommand+=" -c "+str(columnid)+" " gseacommand+=" -p "+str(args.permutation) gseacommand+=" -g \""+args.gmt_file+"\" " diff --git a/tests/test_smoke.py b/tests/test_smoke.py index cdc1cca..15fbc54 100644 --- a/tests/test_smoke.py +++ b/tests/test_smoke.py @@ -225,6 +225,50 @@ def test_gsea_degenerate_pathway_is_not_maximally_significant(tmp_path): assert float(row[4]) == 1.0 +def _gsea_es(tmp_path, rank_text, header_flag): + """Run mageckGSEA on ``rank_text`` and return the enrichment score of PW.""" + rank = tmp_path / f"rank_{'H' if header_flag else 'plain'}.txt" + rank.write_text(rank_text) + gmt = tmp_path / "p.gmt" + # G14..G18 hold the lowest scores, so they cluster at one end of the ranking + # and produce a robustly non-zero enrichment score (~0.87) rather than a + # near-zero value whose textual form differs across platforms. + gmt.write_text("PW\tna\tG14\tG15\tG16\tG17\tG18\n") + out = tmp_path / f"out_{'H' if header_flag else 'plain'}.txt" + cmd = ["mageckGSEA", "-c", "1", "-p", "100", + "-g", str(gmt), "-r", str(rank), "-o", str(out)] + if header_flag: + cmd.append("-H") + result = subprocess.run(cmd, capture_output=True, text=True) + assert result.returncode == 0, result.stderr + rows = { + line.split("\t")[0]: line.split("\t")[2] + for line in out.read_text().splitlines()[1:] + } + return rows["PW"] + + +def test_gsea_skip_header_matches_headerless(tmp_path): + """Regression: ``-H`` must drop the rank-file header, not rank it as a gene. + + ``mageck2 pathway --method gsea`` passes a ``*.gene_summary.txt`` whose first + row is a header. Without skipping it, mageckGSEA parsed the header as a + phantom gene named ``id`` (score ``atof("neg|score") == 0``), inflating the + gene count and shifting every pathway's ranking, percentile, and permutation + math. A headered file read with ``-H`` must therefore score identically to + the same data with no header line at all. + """ + genes = "".join(f"G{i}\t{20 - i}\n" for i in range(1, 21)) # G1..G20, distinct scores + headerless = _gsea_es(tmp_path, genes, header_flag=False) + headered = _gsea_es(tmp_path, "id\tscore\n" + genes, header_flag=True) + + # -H makes the headered input score identically to the headerless data. + assert headered == headerless + # Sanity check that a real (non-degenerate) enrichment score was computed, + # so the equality above is meaningful and not a pair of zeros. + assert float(headered) > 0.5 + + def test_explicit_trim5_does_not_crash(tmp_path): """Regression test for mageck2-doc issue #1.