arnesmits / DEP

DEP package
27 stars 12 forks source link

Adjusted P-values lower than raw P values #21

Closed ScienceBear closed 4 years ago

ScienceBear commented 4 years ago

Hi Arne,

Thanks for the great solution for proteome analysis!

I'm finding that some of my more significant hits are showing lower P-values after adjustment (e.g. 1E-06 to 4E-13. The model I'm fitting is a 3 contrast model.

I dug into DEP's code and it seems it uses limma's default method for adjustment (BH). I've never had this issue with limma before, does DEP take the adjusted P-values from limma and further modify them somehow?

Best wishes,

Luke Dabin

ericloud commented 4 years ago

Hi,

I'm not an author of DEP, but I have your answer. p-value adj. are computed using fdrtool: http://www.strimmerlab.org/software/fdrtool/

I have also some trouble with this method. (see #7) My temporary solution is to modify DEP to disable fdrtool and return p-values adj computed by Limma.

Eric.

ScienceBear commented 4 years ago

Thanks Eric for such a quick response! I expected there was a post-limma step going on. I will take the unadjusted P-values and adjust them using the BH method myself.

m-pauper commented 2 years ago

Hello,

I have also noticed in my case p-adj being lower than original p-values.

Could someone please confirm that fdrtool is not to be trusted? How exactly should the code be modified to use another method?

Thanks!

anissaguillemin commented 1 year ago

Hello! I am also interested to know how exactly the code should be modified. Thanks!

ericloud commented 1 year ago

You don't need to modify DEP, just run DEP and recompute padj column with the method of your choice. Using BH method, it's something like:

df_dep_results$qval = p.adjust(df_dep_results$pvalues, method="BH")

Otherwise to desactivate fdrtool it's here: https://github.com/arnesmits/DEP/blob/b425d8d0db67b15df4b8bcf87729ef0bf5800256/R/functions.R#L975 with something like:

res$qval <- res$adj.P.Val
#fdr_res <- fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE)
#res$qval <- fdr_res$qval
#res$lfdr <- fdr_res$lfdr
anissaguillemin commented 1 year ago

Thanks, I was about to do the second option, but the first one is indeed better.