Skip gene_summary header in mageckGSEA (fixes phantom "id" gene)

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
 

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<string, vector<string> > 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: "<<e.error()<<" for arg "<<e.argId();
@@ -154,6 +158,12 @@ int parserankfile(string rankf,vector<string>& gname, vector<double>& 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++;

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+"\" "

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.