Skip to content

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

Open
ronstewart wants to merge 4 commits into
mainfrom
fix-gprofiler-bugs
Open

Fix rrvgo bugs in gprofiler.r; add multi-ontology reduction, treemaps,#1
ronstewart wants to merge 4 commits 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
ronstewart and others added 3 commits July 27, 2026 16:22
…atting; fix Docker build

- Only shorten term names to their last 3 words when GOBP_-prefixed
  (MSigDB-style); other sources were getting mangled into unreadable
  fragments. Restores this with the correct scoping, replacing the
  blanket disable from the previous commit.
- Diverging/combined comparison output filenames now include both
  dataset names instead of just file1's.
- The skip-clustering fallback for singleton-term ontologies now looks up
  GO.db's space-separated term description instead of reusing the raw,
  underscored source term name, so labels stay consistent with the rest
  of the plot (e.g. "heterotrimeric G-protein complex" instead of
  "heterotrimeric_G-protein_complex").
- Add overlapping apoptosis-related genes to the two gprofiler example
  files so the diverging comparison has a shared, testable term between
  them (mirrors the same test-data change in scRNAseq_downstream).
- Fix deseq2/build/Dockerfile: it was never updated for the rrvgo work in
  the previous commit, so the built image couldn't run this script at all
  (missing rrvgo/GO.db/org.Dr.eg.db/org.Ss.eg.db/devEMF/ggrepel/tidyr/
  tibble/jsonlite/reshape2/dplyr/ggplot2). Also adds libglpk-dev
  (rrvgo's igraph dependency) and libmagick++-dev (GSVA's
  SpatialExperiment dependency), both needed at the system level, and
  drops 'GSVABase', a package name that has never existed on Bioconductor
  and was silently failing on every build.

Verified by rebuilding the Docker image from scratch and running
gprofiler.r end-to-end with --file2 (diverging comparison), reviewing
the generated CSVs/PDFs.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
sign(opt$lfc) errors with "non-numeric argument to mathematical function"
when opt$lfc is NULL, so passing --padj alone (filter by significance
only, no fold-change threshold -- a legitimate combination) crashed the
script. Add the missing filtered-by-padj-only branch, matching
scRNAseq_downstream's gprofiler.r, which already handles this case.

Verified via a standalone logic test (parses cleanly, correctly filters
to the significant row when only --padj is set); didn't need a full
rebuild since this is pure script logic with no new dependencies.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…d output dir

- Add individual raw + rrvgo-reduced bar plots (_barplot.pdf/_reduced_barplot.pdf)
  per input file, matching scRNAseq_downstream's gprofiler.r. Previously the
  only barplots produced were the combined/diverging ones (always named
  "_combined_barplot*"), even when --file2 wasn't given -- there was no
  standalone "just this one file" plot output at all.
- Fix --genes_only/--lower: both are declared type="character" in optparse
  but were checked with `opt$x == TRUE`, which in R only matches the exact
  string "TRUE" -- "true", "T", "1" etc. silently evaluate to FALSE with no
  warning (verified directly). Now uses isTRUE(as.logical(opt$x)), matching
  scRNAseq_downstream's gprofiler.r, which never had this problem.
- Every run now writes into its own output_gprofiler_<timestamp>
  subdirectory under --output_dir instead of directly into --output_dir --
  previously a rerun silently overwrote the previous run's files under the
  same names. Implemented by reassigning opt$output_dir itself right after
  parsing, so every existing write site (including sessionInfo.txt) picks
  up the timestamped path with no other changes needed.
- Rename deseq2/example_data/gprofiler/immune_response_degenes.txt to
  immune_response_markers_degenes.txt (it's a curated marker-gene list,
  not itself a "vs" contrast like its sibling file) and change the
  diverging-comparison combined output's filename join from "_vs_" to
  "_AND_", since --genes' own name can already be a "vs" contrast (e.g.
  photoreceptor_vs_other) -- "_vs_" there read as a confusing nested
  three-way comparison. Mirrors the same rename/fix in
  scRNAseq_downstream.

Verified by rebuilding the local test image and rerunning end-to-end
multiple times: confirmed distinct timestamped output dirs across
consecutive runs, correct per-file barplot content, and correct combined
output naming.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
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