12801 bug in gsea module related with multiple gene sets processing - #12821
12801 bug in gsea module related with multiple gene sets processing#12821brovolia wants to merge 2 commits into
Conversation
| def rpt_label = prefix.replaceAll('\\.$', '') // Remove any trailing dots from prefix when passed as report label, so GSEA doesn't produce double-dotted top-level outputs | ||
| def chip_command = chip ? "-chip $chip -collapse true" : '' | ||
| // gsea-cli crashes if given more than one -gmx file; callers must invoke this process once per gene set file. | ||
| def gmx = gene_sets instanceof List ? gene_sets[0] : gene_sets |
There was a problem hiding this comment.
I don't think silently selecting the first gene set is the right workaround, I think that's worse than erroring when someone supplied the wrong input shape.
We should either explicitly error when a list is supplied, or just remove this and let the fail continue (since any error is at the workflow level)
There was a problem hiding this comment.
Agree, I will drop the instanceof List entirely
There was a problem hiding this comment.
Actually, I think maybe we just need to make a comma-separated version?
| def gmx = gene_sets instanceof List ? gene_sets[0] : gene_sets | |
| def gmx = gene_sets instanceof List ? gene_sets.join(',') : gene_sets |
There was a problem hiding this comment.
No, in this case, gsea-cli will generate an analysis name by appending the list of collection names. I think it should be fixed at the pipeline-level
| def prefix = task.ext.prefix ?: "${meta.id}" | ||
| def rpt_label = prefix.replaceAll('\\.$', '') // Remove any trailing dots from prefix when passed as report label, so GSEA doesn't produce double-dotted top-level outputs | ||
| def chip_command = chip ? "-chip $chip -collapse true" : '' | ||
| // gsea-cli crashes if given more than one -gmx file; callers must invoke this process once per gene set file. |
There was a problem hiding this comment.
We should probably quote the file in the chip_command too
| # Un-timestamp the outputs for path consistency | ||
| mv ${rpt_label}.Gsea.*/* . | ||
| # Only rmdir actual directories: the .rpt file's name also matches "*.Gsea.*". | ||
| mv *.Gsea.*/* . |
There was a problem hiding this comment.
Is it necessary to have made this glob broader?
There was a problem hiding this comment.
yeah.. I will revert it back
pinin4fjords
left a comment
There was a problem hiding this comment.
Ran this for real before going further — pulled the gsea-cli 4.3.2 container and ran the actual script against real GCT/CLS/CHIP data plus two split Hallmark GMT files.
- Comma-joined
-gmx file1.gmt,file2.gmt: succeeds, combines both files (Got gsets: 50), output dir is plain<rpt_label>.Gsea.<timestamp>— no bracketed name list. - Space-separated (what this module currently produces): also exits 0, but gsea-cli silently keeps only the first file and drops the rest (
to parse>hallmark_part1.gmt< got: [hallmark_part1.gmt]). No crash.
So the real bug is silent data loss, not a crash, and join(',') fixes it. Couldn't reproduce the crash or the bracketed-name naming issue — details inline.
| def rpt_label = prefix.replaceAll('\\.$', '') // Remove any trailing dots from prefix when passed as report label, so GSEA doesn't produce double-dotted top-level outputs | ||
| def chip_command = chip ? "-chip $chip -collapse true" : '' | ||
| // gsea-cli crashes if given more than one -gmx file; callers must invoke this process once per gene set file. | ||
| def gmx = gene_sets instanceof List ? gene_sets[0] : gene_sets |
There was a problem hiding this comment.
gsea-cli doesn't crash on multiple files — it silently keeps only the first and drops the rest (confirmed with a real run). It does support multiple files via one comma-joined value (confirmed against GeneSetMatrixMultiChooserParam and a real two-file run combining both, Got gsets: 50).
| def gmx = gene_sets instanceof List ? gene_sets[0] : gene_sets | |
| def gmx = gene_sets instanceof List ? gene_sets.join(',') : gene_sets |
There was a problem hiding this comment.
ok, so here is the question then, do you want one GSEA_GSEA call per contrast with multiple slurm jobs or a multi-collection call per contrast so gsea-cli will handle it? If you want the second option, then I guess conf/modules.config should also be corrected after
There was a problem hiding this comment.
I'm not sure there's a good reason to separate them into multiple jobs?
But in any case, it's a slightly orthogonal question. From the module's point of view as a component it should be able to deal with the multiple-input method, since the GSEA CLI can.
| def prefix = task.ext.prefix ?: "${meta.id}" | ||
| def rpt_label = prefix.replaceAll('\\.$', '') // Remove any trailing dots from prefix when passed as report label, so GSEA doesn't produce double-dotted top-level outputs | ||
| def chip_command = chip ? "-chip $chip -collapse true" : '' | ||
| // gsea-cli crashes if given more than one -gmx file; callers must invoke this process once per gene set file. |
There was a problem hiding this comment.
This comment's wrong once gmx is comma-joined — worth dropping.
| @@ -45,26 +45,32 @@ process GSEA_GSEA { | |||
| def prefix = task.ext.prefix ?: "${meta.id}" | |||
| def rpt_label = prefix.replaceAll('\\.$', '') // Remove any trailing dots from prefix when passed as report label, so GSEA doesn't produce double-dotted top-level outputs | |||
| def chip_command = chip ? "-chip $chip -collapse true" : '' | |||
There was a problem hiding this comment.
Since everything else is getting quoted:
| def chip_command = chip ? "-chip $chip -collapse true" : '' | |
| def chip_command = chip ? "-chip \"$chip\" -collapse true" : '' |
| # Only rmdir actual directories: the .rpt file's name also matches "*.Gsea.*". | ||
| mv *.Gsea.*/* . | ||
| for gsea_dir in *.Gsea.*/; do | ||
| [ -d "\$gsea_dir" ] && rmdir "\$gsea_dir" | ||
| done |
There was a problem hiding this comment.
Tried to reproduce this naming problem — a real comma-joined -gmx across two GMT files produced a plain test_comma3.Gsea.<timestamp> dir, which the original glob already matches. Checked ToolReport.generateReportDir in gsea-desktop v4.3.2 too — built from --rpt_label + tool name + timestamp only. Suggest reverting once gmx is comma-joined:
| # Only rmdir actual directories: the .rpt file's name also matches "*.Gsea.*". | |
| mv *.Gsea.*/* . | |
| for gsea_dir in *.Gsea.*/; do | |
| [ -d "\$gsea_dir" ] && rmdir "\$gsea_dir" | |
| done | |
| mv ${rpt_label}.Gsea.*/* . |
I narrowed the removal pattern and used quotation marks, so the white space will not break the processes