Using the celltype enrichment results gen_overlap.symptoms.rds, add an arg to MultiEWCE::prioritise_targets to let users filter by significance.
https://github.com/neurogenomics/MultiEWCE/releases/tag/v0.0.1

E.g. to enable you to make neurogenomics/rare_disease_celltyping#2

Only occurs on windows:
https://github.com/neurogenomics/MultiEWCE/actions/runs/6798449226/job/18482623120#step:2:2571

Run options(crayon.enabled = TRUE, timeout=Sys.getenv("timeout"))
── R CMD build ─────────────────────────────────────────────────────────────────
* checking for file 'D:\a\MultiEWCE\MultiEWCE/DESCRIPTION' ... OK
* preparing 'MultiEWCE':
* checking DESCRIPTION meta-information ... OK
* installing the package to build vignettes
* creating vignettes ... ERROR
Error: --- re-building 'MultiEWCE.Rmd' using rmarkdown
Analysing: HP:0025114 (1/3): 20 genes.
Analysing: HP:0002948 (2/3): 70 genes.
Analysing: HP:0410219 (3/3): 6 genes.
The magick package is required to crop "C:/Users/runneradmin/AppData/Local/Temp/Rtmpwf8x0G/Rbuild11c444b82d73/MultiEWCE/vignettes/MultiEWCE_files/figure-html/display-1.png" but not available.
Quitting from lines 135-136 [merge] (MultiEWCE.Rmd)
Error: Error: processing vignette 'MultiEWCE.Rmd' failed with diagnostics:

In MultiEWCE 0.1.8 i made a number of upgrades/fixes:

New features

• MultiEWCE finally gets a hex sticker!
• gen_results / gen_overlap
  • Check for existing results and import if already there.
  • Name all results "gen_results.rds" or "gen_overlap.rds" to avoid
    rerunning duplicate analyses on HPC.
• Update rworkflows.yml

Bug fixes

• Fix unit tests and examples to use "hpo_id" instead of "hpo_name".
• load_hpo_graph:
  • export
  • Regenerate and update "hpo_graph.rds" file.
• Drastically reduce time to run examples.
• ontology_plot
  • Fix function and add test.
• DESCRIPTION
  • Depends: R (>= 2.10) --> Depends: R (>= 4.1),
    to ensure |> function available.
  • Rewrite Description field to reflect MultiEWCE's current purpose.

Remaining bugs

Everything seems to pass fine locally, but getting a couple of test failures on GHA:
https://github.com/neurogenomics/MultiEWCE/actions/runs/6755745935/job/18364474266

May have to do with file version updates in HPOExplorer and some caching happening on my local machine.

Saving plot ==> /tmp/RtmpEDtTTF/working_dir/RtmpQIWJbO/filec13b3431c62_prioritise_targets_network.html
[ FAIL 4 | WARN 2 | SKIP 0 | PASS 24 ]

Error: R CMD check found ERRORs
Execution halted
══ Failed tests ════════════════════════════════════════════════════════════════
── Failure ('test-ewce_para.R:45:3'): ewce_para works ──────────────────────────
sum(list_names %in% unique(all_results1$hpo_name)) is not more than length(list_names) - 1. Difference: -1
── Failure ('test-ewce_para.R:47:3'): ewce_para works ──────────────────────────
sum(list_names %in% unique(all_results2$hpo_name)) is not more than length(list_names) - 1. Difference: -1
── Failure ('test-gen_results.R:13:3'): gen_results works ──────────────────────
sum(list_names %in% unique(all_results$hpo_id)) is not more than length(list_names) - 1. Difference: -1
── Failure ('test-gen_results.R:15:3'): gen_results works ──────────────────────
nrow(all_results[q <= 0.05, ]) is not more than 10. Difference: -2

[ FAIL 4 | WARN 2 | SKIP 0 | PASS 24 ]
Error: Test failures
Execution halted

Session info

</details>

@NathanSkene @KittyMurphy
I just noticed that the CTD originally included within MultiEWCE is incorrect. It was not processed to the level of gene properly.

MultiEWCE::load_example_ctd()

Instead, it appears to be at the level of transcripts, which were then randomly assigned unique gene-level IDs. You can also tell this right away when looking at the number of dimensions in the CTD.

@bobGSmith can you please confirm which version of the CTD you used in your original EWCE analyses? If it was this version, I will need to redo all of the analyses. Otherwise, our gene counts will be underestimated by ~2/3 and thus the celltype signatures would be degraded.

An important question for any enrichment-based analysis is: how do we determine the appropriate set of background genes?

In the case of rare disease celltyping project, I think there's two potential backgrounds that could be used.

Options

1. All HPO genes

All genes included in the HPO gene annotations.
This is essentially saying "is this phenotype's genes enriched in a celltype relative to rare disease genes in general?"

> gene_data <- HPOExplorer::load_phenotype_to_genes()
Reading cached RDS file: phenotype_to_genes.txt
+ Version: v2023-10-09
> length(unique(gene_data$gene_symbol))
[1] 5005

2. All human genes

All genes in the human genome, or at least those that appear in the CTD. When the CTD is from mouse, this background would be further reduced to those that have human 1:1 orthologs.

This is essentially saying "is this phenotype's genes enriched in a celltype relative to all genes that are expressed and informative in some cell type?"

genelistSpecies = "human"
sctSpecies = "human"
> bg = orthogene::create_background(
+     species1 = genelistSpecies,
+     species2 = sctSpecies)
Retrieving all genes using: homologene.
Retrieving all organisms available in homologene.
Mapping species name: human
Common name mapping found for human
1 organism identified from search: 9606
Gene table with 19,129 rows retrieved.
Returning all 19,129 genes from human.
Returning 19,129 unique genes from entire human genome.

> length(bg)
[1] 19129

Internally, EWCE will use the CTD genes to subset this list further, as it can only sample specificity scores from genes that are included in the CTD:

• https://github.com/NathanSkene/EWCE/blob/f2d0f49eae1718862a3abfdaed73db00c0b93c18/R/bootstrap_enrichment_test.R#L136
• https://github.com/NathanSkene/EWCE/blob/master/R/get_summed_proportions.r

Conclusions

@NathanSkene @KittyMurphy I think we concluded the latter is the way to go, especially as this is concordant with how EWCE behaves by default (unless the user supplies their own custom bg= arg).

Just making a note as I figured out how to solve this but haven't had time to change it on the web app. You need to create this function

phenotype_namer <- function(HPO_id) {
return (hpo$name[HPO_id])
} 

Then when you call the onto_plot function, assign the phenotype_name function to the label argument e.g.

onto_plot(hpo,  
terms= terms ,
label=phenotype_namer,  
shape="circle", 
edge_attributes = list(color="black"),  
fontsize = 80)

Currently I'm using a table from Momoko's thesis to filter the phenotypes that are in highest severity tier. But can try basing this off of the annotation data provided by HPO.

Parallelising over groups in R (e.g. each phenotype's gene list) can use a ton of memory, because you can accidentally copy all objects used inside the loop, which multiplies your memory usage when working with large datasets.

gene_data <- HPOExplorer::load_phenotype_to_genes()
gene_data[,HPO_ID.LinkID:=paste(HPO_ID,LinkID,sep=".")]
nrow(gene_data) # 968,116
length(unique(gene_data$HPO_ID.LinkID)) #621,906

Even on large-memory machines like the Threadripper (250Gb) this quickly gets overloaded bc each parallel process uses 17.6Gb of memory to store the entire gene_data obj over and over again.
17.6 * 50 cores = 835Gb of memory! This means that parallelising chokes up the threadripper's memory and slows it down to a grinding halt, making processing the data even slower than if you had just single-threaded it. This is exactly what happened when i tried this version of the gen_overlap func:
https://github.com/neurogenomics/MultiEWCE/blob/25d26a41096902607a4f595343f2f585dad9f819/R/gen_overlap.R#L49

A better way might be to split the data first into chunks, and then iterate over the chunks:

 split.data.table <- utils::getFromNamespace("split.data.table","data.table")
gene_data_split <- split.data.table(x = gene_data,
                                    by ="HPO_ID.LinkID")

Related posts

https://stackoverflow.com/questions/19082794/speed-up-data-table-group-by-using-multiple-cores-and-parallel-programming
https://stackoverflow.com/questions/14759905/data-table-and-parallel-computing

👇 Contains some useful benchmarking
Rdatatable/data.table#2223

Killing zombie processes

This happens when you launch some paralleised func in R and then decide to stop it midway. A bunch of "zombie" processes are leftover.

Restarting the R session

Sometimes this works, other times not so much. Might not work when memory is so overloaded that you can't restart the R session.

via htop

Didn't seem to work.

Via inline

Someone suggested this, but didn't seem to do anything for me on the Threadripper.

https://stackoverflow.com/questions/25388139/r-parallel-computing-and-zombie-processes

via future

Using future instead of parallel might give me more control over this, but it has to done before launching the processes.
HenrikBengtsson/future#93

via Docker

To kill them, the only effective way I've found is to restart the container:

docker stop <id>
docker start <id>

@Al-Murphy probably relevant to you too

prioritise_targets_network uses visNetwork to generate really nice interactive plots like this:

Unfortunately, visNetwork doesn't provide any means of exporting the html/javascript-rendered networks as PDFs or PNGs at sufficiently high resolution. The "export PDF" feature is basically just saving a screenshot...
A version of this feature that actually works has been requested many times, but apparently it's a bit tricky and the developer hasn't had time to address it:
https://github.com/datastorm-open/visNetwork/issues?q=visExport

My options include:

• Hacking a solution to export a decent PNG/PDF for publication figures.
• Rewriting the entire network plot code so that it uses something else, like ggnetwork. But it will take some time for me to figure out how to make this look decent.

This now results a named list, one of which is "report":

res <- MultiEWCE::prioritise_targets()

"report" is a recording of how many phenotypes/celltypes/genes remained after each filtering step.
@NathanSkene This will help us to figure out where phenotypes are getting dropped the most.

I've tried expanding the Tier assignments using a combination of Momoko's curated + my automated approaches (described here #5), but this still seems to be a step where the most phenotypes are dropped. Not necessarily a problem, just something to be aware of.

Seems to only be for certain cell types, eg.

Reprex

res=MultiEWCE::ggnetwork_plot_full( 
    results = load_example_results(),  
    cell_type = "Chromaffin cells", 
    q_threshold = 0.005, 
    fold_threshold = 1) 

Output

ggnetwork_plot_full
Subsetting results by phenotype.
get_cell_ontology
Importing existing file: /Users/schilder/Library/Caches/org.R-project.R/R/HPOExplorer/data/phenotype_to_genes.txt
5 associations remain after filtering.
Creating adjacency matrix for 5 terms.
Making hoverboxes from: 'Phenotype', 'HPO_ID', 'description', 'CellType', 'p', 'q', 'fold_change'
Making phenotype network object.
Getting ontology level for 5 HPO IDs.
Creating ggnetwork plot.
Error in cbind(x, y, others) : object 'x' not found

Full traceback

Error in cbind(x, y, others) : object 'x' not found
14.
cbind(x, y, others)
13.
rbind(cbind(x, y, others), cbind(x = xend, y = yend, others))
12.
eval(substitute(expr), data, enclos = parent.frame())
11.
eval(substitute(expr), data, enclos = parent.frame())
10.
with.default(data, {
rbind(cbind(x, y, others), cbind(x = xend, y = yend, others))
})
9.
with(data, {
rbind(cbind(x, y, others), cbind(x = xend, y = yend, others))
})
8.
to_basic.GeomSegment(data[[i]], prestats_data[[i]], layout, params[[i]],
p)
7.
to_basic(data[[i]], prestats_data[[i]], layout, params[[i]],
p)
6.
layers2traces(data, prestats_data, layout, plot)
5.
gg2list(p, width = width, height = height, tooltip = tooltip,
dynamicTicks = dynamicTicks, layerData = layerData, originalData = originalData,
source = source, ...)
4.
ggplotly.ggplot(p = network_plot, tooltip = tooltip)
3.
plotly::ggplotly(p = network_plot, tooltip = tooltip) at ggnetwork_plot.R#61
2.
plotly::layout(plotly::ggplotly(p = network_plot, tooltip = tooltip),
hoverlabel = list(align = "left")) at ggnetwork_plot.R#61
1.
HPOExplorer::ggnetwork_plot(phenoNet = phenoNet, colour_var = colour_var,
interactive = interactive, verbose = verbose)

Session info

</details>

Count the number of unique disease IDs each phenotype in the HPO is associated with, and store this info in a new column called "n_diseases":
https://github.com/neurogenomics/HPOExplorer/blob/master/R/add_ndisease.R

I now use this information to filter phenotypes associated with >20 diseases (by default). User can adjust this with the new arg: prioritise_targets(pheno_ndiseases_threshold=)

ie. only label celltype nodes with celltype variables, only label phenotypes with HPO phenotype variables, etc.

In prioritise_targets, include an option to only keep celltype-phenotype associations that appear in the matching celltype-symptom (phenotype + disease) results.

This means that there only genes identified as putative mediator genes in celltype-symptom tests will be kept. Any phenotypes that don't have any genes left after this filter will be dropped.

More granular control can be achieved using the symptom_p_threshold arg as well.
#9