I am trying to reproduce a functional enrichment analysis with Mummichog from the web app with the R package. For this I have to set a p-value cutoff with the `SetMummichogPval()` function. As standard settings, the web app claims to use the top 10 % peaks. To reproduce this in R, I ranked the features according to their p-values in an ascending order (lowest p-values/highest significance at the top) and I took the lowest value from the most significant 10 % (= the 10th percentile) as the cutoff. When I did it like this, the value I got was around 0.008. The web app, which claims to do the same thing, gives me a cutoff of 0.005 and with this value I get slightly different results from the enrichment analysis. I concluded that the value is somehow rounded in the web app in a way that is not documented (at least I couldn’t find it). I played around with different input files but could not figure out exactly how the rounding is done.
Could you please explain how the rounding is done in the web app?
You are comparing a single step out of the analysis context, please refer to ‘R History’ to see the exact steps used in the web-based version.
Thank you very much for your answer!
Unfortunately, the Rhistory does not state how the cutoff is derived. It is just set to that value without further explanation. Below is the full Rhistory. The marked value is the parameter in question.
mSet<-InitDataObjects(“mass_all”, “mummichog”, FALSE, 150)
mSet<-SetPeakFormat(mSet, “rmp”)
mSet<-UpdateInstrumentParameters(mSet, 5.0, “mixed”, “yes”, 0.02);
mSet<-Read.PeakListData(mSet, “Replacing_with_your_file_path”);
mSet<-SanityCheckMummichogData(mSet)
mSet<-SetPeakEnrichMethod(mSet, “mum”, “v2”)
mSet<-SetMummichogPval(mSet, 0.005)
mSet<-PerformPSEA(mSet, “hsa_kegg”, “current”, 3 , 100)
mSet<-PlotPeaks2Paths(mSet, “peaks_to_paths_0_”, “png”, 150, width=NA)
The interface of the web app claims that it corresponds to the top 10% of peaks but it is not the case. If it was the case, the cutoff would be around 0.008. In my case, this means that out of approximately 6500 peaks, only around 450 are considered significant and not 650 which I would expect. This leads to different results when I run the enrichment analysis.
Again, thank you for your time!
R command history is mainyl for transparency and reproducibility, as web interface can not show everything. If you are familiar with R, you can view the underlying code and method documentation.