Skip to content

Fix rrvgo bugs in gprofiler.r; add multi-ontology reduction, treemaps, - #1

Open
ronstewart wants to merge 1 commit into
mainfrom
fix-gprofiler-bugs
Open

Fix rrvgo bugs in gprofiler.r; add multi-ontology reduction, treemaps,#1
ronstewart wants to merge 1 commit into
mainfrom
fix-gprofiler-bugs

Conversation

@ronstewart

Copy link
Copy Markdown

JSON hierarchy, and significance-colored diverging plots

Ported over the bug fixes and later feature work done on scRNAseq_downstream's copy of this script (which started from this one), adapted to this script's CLI-args + optional two-file diverging-comparison design (--file2/--color1/--label1).

Three bug fixes (independent of everything else below)

  • Log-base inconsistency: the raw bar plot used natural log (-log(p.adj)) while the rrvgo score used log10 (-log10(p.adj)), so the same term's bar height looked wildly different between the two plots (~15 vs ~6.5 for the same p-value) despite the axis label already claiming log10 for both. Both now use log10 consistently.
  • calculateSimMatrix()'s last internal step is m[!out, !out], which R's default matrix indexing silently collapses to a bare scalar whenever exactly one term survives orgdb/ancestor filtering -- losing the dimnames that would identify which term it was. A scalar NA means truly zero terms mapped; a non-NA scalar means exactly one term mapped. The previous check (!is.null(dim(sim_matrix)) alone) treated both cases as "zero mapped," silently dropping a real, valid single-term cluster.
  • orgdb was hardcoded to "org.Dr.eg.db" regardless of --organism, so e.g. running with --organism hsapiens would still (incorrectly) use zebrafish annotations for the rrvgo step. Now looked up from --organism via deseq2/build/src/organism_orgdb_map.txt (same design as scRNAseq_downstream's data/organism_orgdb_map.txt -- add a row there and install the matching Bioconductor package for a new organism). Resolved relative to this script's own location (not getwd()), since this script has no "run from repo root" convention.

Feature parity with scRNAseq_downstream, adapted for --file2

All four of these live inside run_gprofiler_analysis(), which already runs once per gene-list file -- so file1 and file2 (when given) each get their own full set of per-file outputs independently; only the final combined diverging plots needed new logic to merge them.

  • Multi-ontology support: new --ontologies flag (comma-separated, e.g. "BP" or "BP,MF,CC"; defaults to "BP" to match prior behavior). Semantic similarity is only meaningful within one GO ontology at a time, so each requested ontology is run through calculateSimMatrix/reduceSimMatrix separately and stacked with an ontology column. Terms whose ontology can't be resolved, or whose ontology isn't in --ontologies, or that calculateSimMatrix couldn't map in orgdb, are recorded (with why) in a new _reduced_excluded.csv per file instead of silently vanishing.
  • Significance-colored diverging plots: the combined plots used to color bars by file1-vs-file2 identity (opt$color1 vs red), which meant they couldn't also show significance, and conflicted with wanting a continuous color scale like scRNAseq_downstream's. Since bar position (left/right of zero) already encodes file1 vs file2 -- that's the whole point of the sign flip -- fill color is now a continuous -log10(p.adj) gradient (blue->red) on both sides instead, applied to a new unsigned significance column carried alongside the signed score (used only for position). No information is lost: the y-axis label already spells out which side is R/NR.
  • reduced_barplot is now faceted by ontology (one panel per BP/MF/CC), and grouped by parentTerm rather than individual term, matching scRNAseq_downstream -- a term absorbed into a more significant parent's cluster no longer gets its own bar.
  • Per file, per ontology with data: reduced_treemap.pdf (rrvgo's own cluster visualization -- doesn't facet, so one PDF per ontology instead of trying to force them into one), plus _reduced_grouping.csv (one row per parent cluster, listing member terms directly) and _reduced_hierarchy.json (the same grouping as nested JSON: ontology -> parent cluster -> member terms, browsable in a JSON viewer's fold view -- easier to scan than the flat CSV or a busy treemap with many clusters).

Test data + .gitignore

Added deseq2/example_data/gprofiler/ (mirroring the existing deseq2/example_data/ convention for deseq2.r's own test inputs):

  • photoreceptor_vs_other_degenes.txt: real zebrafish photoreceptor marker genes, converted from scRNAseq_downstream's data/test_gprofiler/photoreceptors_vs_other_lfc0.25.txt (Seurat avg_log2FC/p_val_adj column names -> this script's expected DESeq2-style log2FoldChange/padj).
  • immune_response_degenes.txt: synthetic (clearly marked as such in its header comment) immune/inflammatory-response gene data, specifically to pair with the photoreceptor file as --genes/--file2 for a genuinely asymmetric diverging-plot test -- testing --file2 against the same file twice is trivially symmetric and doesn't exercise the position/color encoding logic above.

Both organism_orgdb_map.txt and the new example_data files were silently excluded by .gitignore's existing blanket *.txt rule (no .txt file has ever been tracked in this repo). Added targeted ! exceptions for these specific paths rather than removing the blanket rule, so future generated/log .txt output stays ignored.

Verified end-to-end: single-file runs (multi-ontology, all new per-file outputs) and two-file diverging runs, including a run against real photoreceptor vs. synthetic immune data confirming genuinely asymmetric results (different terms, different magnitudes, no forced symmetry) with correct position-for-direction/color-for-significance rendering.

JSON hierarchy, and significance-colored diverging plots

Ported over the bug fixes and later feature work done on
scRNAseq_downstream's copy of this script (which started from this one),
adapted to this script's CLI-args + optional two-file diverging-comparison
design (--file2/--color1/--label1).

## Three bug fixes (independent of everything else below)

- Log-base inconsistency: the raw bar plot used natural log
  (-log(p.adj)) while the rrvgo score used log10 (-log10(p.adj)), so the
  same term's bar height looked wildly different between the two plots
  (~15 vs ~6.5 for the same p-value) despite the axis label already
  claiming log10 for both. Both now use log10 consistently.
- calculateSimMatrix()'s last internal step is `m[!out, !out]`, which R's
  default matrix indexing silently collapses to a bare scalar whenever
  exactly one term survives orgdb/ancestor filtering -- losing the
  dimnames that would identify which term it was. A scalar NA means truly
  zero terms mapped; a non-NA scalar means exactly one term mapped. The
  previous check (`!is.null(dim(sim_matrix))` alone) treated both cases
  as "zero mapped," silently dropping a real, valid single-term cluster.
- orgdb was hardcoded to "org.Dr.eg.db" regardless of --organism, so
  e.g. running with --organism hsapiens would still (incorrectly) use
  zebrafish annotations for the rrvgo step. Now looked up from
  --organism via deseq2/build/src/organism_orgdb_map.txt (same design as
  scRNAseq_downstream's data/organism_orgdb_map.txt -- add a row there
  and install the matching Bioconductor package for a new organism).
  Resolved relative to this script's own location (not getwd()), since
  this script has no "run from repo root" convention.

## Feature parity with scRNAseq_downstream, adapted for --file2

All four of these live inside run_gprofiler_analysis(), which already
runs once per gene-list file -- so file1 and file2 (when given) each get
their own full set of per-file outputs independently; only the final
combined diverging plots needed new logic to merge them.

- Multi-ontology support: new --ontologies flag (comma-separated, e.g.
  "BP" or "BP,MF,CC"; defaults to "BP" to match prior behavior).
  Semantic similarity is only meaningful within one GO ontology at a
  time, so each requested ontology is run through
  calculateSimMatrix/reduceSimMatrix separately and stacked with an
  `ontology` column. Terms whose ontology can't be resolved, or whose
  ontology isn't in --ontologies, or that calculateSimMatrix couldn't
  map in orgdb, are recorded (with why) in a new
  <name>_reduced_excluded.csv per file instead of silently vanishing.
- Significance-colored diverging plots: the combined plots used to color
  bars by file1-vs-file2 identity (opt$color1 vs red), which meant they
  couldn't also show significance, and conflicted with wanting a
  continuous color scale like scRNAseq_downstream's. Since bar
  *position* (left/right of zero) already encodes file1 vs file2 --
  that's the whole point of the sign flip -- fill color is now a
  continuous -log10(p.adj) gradient (blue->red) on both sides instead,
  applied to a new unsigned `significance` column carried alongside the
  signed `score` (used only for position). No information is lost: the
  y-axis label already spells out which side is R/NR.
- reduced_barplot is now faceted by ontology (one panel per BP/MF/CC),
  and grouped by parentTerm rather than individual term, matching
  scRNAseq_downstream -- a term absorbed into a more significant
  parent's cluster no longer gets its own bar.
- Per file, per ontology with data: <name>_reduced_treemap_<ontology>.pdf
  (rrvgo's own cluster visualization -- doesn't facet, so one PDF per
  ontology instead of trying to force them into one), plus
  <name>_reduced_grouping.csv (one row per parent cluster, listing member
  terms directly) and <name>_reduced_hierarchy.json (the same grouping as
  nested JSON: ontology -> parent cluster -> member terms, browsable in a
  JSON viewer's fold view -- easier to scan than the flat CSV or a busy
  treemap with many clusters).

## Test data + .gitignore

Added deseq2/example_data/gprofiler/ (mirroring the existing
deseq2/example_data/ convention for deseq2.r's own test inputs):
- photoreceptor_vs_other_degenes.txt: real zebrafish photoreceptor marker
  genes, converted from scRNAseq_downstream's
  data/test_gprofiler/photoreceptors_vs_other_lfc0.25.txt (Seurat
  avg_log2FC/p_val_adj column names -> this script's expected DESeq2-style
  log2FoldChange/padj).
- immune_response_degenes.txt: synthetic (clearly marked as such in its
  header comment) immune/inflammatory-response gene data, specifically to
  pair with the photoreceptor file as --genes/--file2 for a *genuinely*
  asymmetric diverging-plot test -- testing --file2 against the same file
  twice is trivially symmetric and doesn't exercise the position/color
  encoding logic above.

Both organism_orgdb_map.txt and the new example_data files were silently
excluded by .gitignore's existing blanket `*.txt` rule (no .txt file has
ever been tracked in this repo). Added targeted `!` exceptions for these
specific paths rather than removing the blanket rule, so future
generated/log .txt output stays ignored.

Verified end-to-end: single-file runs (multi-ontology, all new per-file
outputs) and two-file diverging runs, including a run against real
photoreceptor vs. synthetic immune data confirming genuinely asymmetric
results (different terms, different magnitudes, no forced symmetry) with
correct position-for-direction/color-for-significance rendering.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
@ronstewart
ronstewart requested a review from bmmoore43 July 24, 2026 16:21
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant