Entering edit mode
5.9 years ago
Biologist
▴
290
I'm doing differential analysis between two groups A & B
With Fold change > 2 and FDR < 0.05 I have got more than 600 differentially expressed genes. Among them I see the following gene.
LogFC Unshrunk.logFC logCPM PValue FDR
RMRP 2.802567464 2.802628518 11.43969321 1.94E-06 2.07E-05
But with Fold change > 5 and FDR < 0.05 I did not find this gene differentially expressed.
Similarly in other analysis between C & D
with Fold change > 2 and FDR < 0.05
LogFC Unshrunk.logFC logCPM PValue FDR
RMRP 9.269439542 9.275584319 12.44069535 1.73E-28 1.48E-26
with Fold change > 5 and FDR < 0.05
LogFC Unshrunk.logFC logCPM PValue FDR
RMRP 9.269439542 9.275584319 12.44069535 1.77E-24 8.20E-22
My question why in the analysis between A&B with Foldchange > 5 I didn't find that gene differentially expressed?
In the analysis between C & D, with Fold change > 2 and Fold Change > 5 I see the "RMRP" gene showing similar values but FDR is changing. Does FDR value change based on Fold change cutoff?
RMRP has a logFC of 2.8 (assuming this is log2) it is 2^2.8 = > 6 fold change. So you are doing something wrong. Did you select logFC > 2?
It is not clear what kind of test you are doing? Is it limma? Or are you using treat?
logFC 2.8 => 6 fold change, which means it is Fold Change > 2 and should also present with Fold Change > 5 also right.
I'm using edgeR, glmTreat function.
glmTreat does not work with a strict cut off, but uses an alternative hypothesis testing. In a normal test the hypothesis is 0 (FC) difference, but with glmTreat you test directly whether there is a 5 fold change. If this is not the case the 5 fold change is not significant (even if you find logFC > 2.8). This logFC is a rough estimate but is not significant via the treat test.
I guess this could be due to multiple testing?
No it is how the test works, read here the details of the test. If you test for > 5 FC, your RMRP gene in your first example doesn't meet the threshold statistically. That the logFC is > 6 FC is not the same, this logFC is an estimate which is not accurate, but more of an indication (they also explain this in the edgeR manual, that logFC are not accurate).
From limma manual
In what context? For edgeR glmTreat? I doubt it. A fold change cut off in limma is something else than glmTreat in edgeR!
If you look at the last two examples in the original question, you can see that uncorrected p-value is also changing, not only the FDR corrected ones.
@OP: I think your question is answered here: https://support.bioconductor.org/p/110243/
To gain a full appreciation of what's going on, you should be outputting, in addition, things like the mean expression, median expression, variance, covariance, and dispersion. I had never assumed that the FDR might be directly 'related to' (correlated with) the magnitude of the fold change. Indeed, the FDR calculation has no direct relationship to the fold-change value - they are calculated independently.
You must also consider the type of normalisation that you've done, and the origin of your counts and how they were 'counted'.