Fix rrvgo bugs in gprofiler.r; add multi-ontology reduction, treemaps, - #1
Open
ronstewart wants to merge 1 commit into
Open
Fix rrvgo bugs in gprofiler.r; add multi-ontology reduction, treemaps,#1ronstewart wants to merge 1 commit into
ronstewart wants to merge 1 commit into
Conversation
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>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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)
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.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.
ontologycolumn. 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.significancecolumn carried alongside the signedscore(used only for position). No information is lost: the y-axis label already spells out which side is R/NR.Test data + .gitignore
Added deseq2/example_data/gprofiler/ (mirroring the existing deseq2/example_data/ convention for deseq2.r's own test inputs):
Both organism_orgdb_map.txt and the new example_data files were silently excluded by .gitignore's existing blanket
*.txtrule (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.