Add calculation of gene set metrics module - #24
Merged
Conversation
|
sjspielman
reviewed
Sep 2, 2026
sjspielman
left a comment
sjspielman
left a comment
There was a problem hiding this comment.
Looks good to me overall, just a few comments with one design suggestion to allow the adjusted p-value to be specified here.
Also a note that I used this script for reviewing: https://github.com/AlexsLemonade/ews-nf/blob/main/modules/metaprograms/resources/usr/bin/04-metaprogram-metrics.R (the links in the PR went to this notebook)
|
|
||
| # run ORA using the top genes | ||
| ora_results <- mp_top_list |> | ||
| purrr::map(function(genes) { |
There was a problem hiding this comment.
Probably we want \(x) {} syntax here
Suggested change
| purrr::map(function(genes) { | |
| purrr::map(\(genes) { |
Member
Author
There was a problem hiding this comment.
Just a note that we can't do this with Nextflow templates. They complain about this syntax which is why we use function here.
| # combine ora results into a single dataframe | ||
| # only include genesets with pvalue < 0.05, the default by enricher | ||
| ora_results_df <- ora_results |> | ||
| purrr::map(function(ora) ora@result) |> |
There was a problem hiding this comment.
Suggested change
| purrr::map(function(ora) ora@result) |> | |
| purrr::map(\(ora) ora@result) |> |
Comment on lines
+43
to
+44
| # returns a df with all gene sets for each metaprogram with p.adjust <= 0.05 | ||
| run_ora <- function(mp_top_list, gene_universe, term2gene_df) { |
There was a problem hiding this comment.
A thought I just had...It might be good to make the adjusted p-value an argument to scripts in this module overall (and eventually, not now, a workflow parameter) in case someone wants to vary for the own usage this e.g. to 0.01
| # keep the total number of genesets the same | ||
| shuffled_ora_results <- mp_geneset_counts[["geneset_count"]] |> | ||
| purrr::set_names(mp_geneset_counts[["metaprogram"]]) |> | ||
| purrr::map(function(count) { |
There was a problem hiding this comment.
Suggested change
| purrr::map(function(count) { | |
| purrr::map(\(count) { |
|
|
||
| # get the top X number of genes for each MP | ||
| mp_top_list <- mp_list |> | ||
| purrr::map(function(weights){ |
There was a problem hiding this comment.
Suggested change
| purrr::map(function(weights){ | |
| purrr::map(\(weights){ |
Co-authored-by: Stephanie J. Spielman <stephanie.spielman@gmail.com>
Member
Author
|
Okay this is now working correctly. I guess pixi doesn't like to build certain bioconductor packages so we have to just use the |
sjspielman
approved these changes
Sep 2, 2026
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
2 participants
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.
This PR adds the second half of the module for calculating metrics, which includes metrics specific to genesets. I took any code related to gene set metrics, including running ORA, out of the original script and added it to the template here.
Similar to the other metrics module, I decided to output TSV files here rather than an RDS since I think it will be easier to work with. I saved the TSV with the full ORA results, the metrics for each metaprogam (gene set specificity), and the background stats.
The other main difference here is that I read in a pre-prepared term2gene file with all genes in the gene set rather than pulling from
msigdbevery time the script is run.The original script also included extracting the top 200 genes from the metaprograms. I think it makes more sense to do this once and store it in the original RDS object that's saved as output from
generate-metaprograms. So I updated that script to include the top genes.When reviewing, please compare implementation to the previous script and look for any discrepancies that might not be intentional as described above.
Note that I'm still testing this since I'm getting an error that there's no package
GO.db. I'm working on fixing this with the environment file, but figured it was still worth getting review started on the rest of it.