Hello there,
I am a little confused by the toptable output (well really the FDR) and thought I'd quickly check my interpretation.
For the following code:
topall <- toptable(fit, coef="CaseAD", genelist=fit$genes, number=Inf, adjust="fdr", p.value=0.01, lfc=log2(1))
Would the interpretation be:
Using a threshold of 0.05 for FDR adj. p.values and a minimal fold change of 2, this analysis produces n number of genes?
And does the above mean that I would look at the adjusted p.values?
I am interested in following MAQC guidelines for this analysis, so would rather restrict my results to a FDR = 0.05 and p.values < 0.01. However, I am not entirely sure, how to get that, since I am not sure what the "adjust=fdr" relates to here?
Any advice?
Many thanks!!!
Thank you!
Just to clarify: Does that mean my lfc is redundant?
Does that also mean you should plot adjusted p.values instead of p.values?
Yes, your log-fold change threshold is redundant in your original post.
For plotting purposes, it's better to use the
-log10(p.value)
, as the enforced monotonicity of the BH correction leads to stretches of genes with the same FDR value, which is not aesthetically pleasing. (For plotting only; calling of DE genes is still done with the FDR.) If you want to visualize the threshold of significance, simply use the largest p-value corresponding to a FDR that is below your threshold. Also check outvolcanoplot
if that's what you intend to make.Thanks Aaron.
Just to clarify, your lfc is redundant because you set lfc = log2(1) = 0, which is the default and imposes no fold change thresholds on your results, so you are only selecting based on adj. p-values. A fold change of 2 as you indicated would be lfc = 1, but see warning in ?topTable on setting p-val and lfc cutoffs simultaneously, see ?treat and ?topTreat for fold change thresholding.