I have an RNA-seq experiment where I have 6 control and 6 case samples. I am getting incredibly small p-value for my DEseq, diferential expression analysis (p-value < 1e-139 for one detector).
nbinomTest(cds[c("xxx1","xxx2")],"Control", "Case")
id baseMean baseMeanA baseMeanB foldChange log2FoldChange
1 xxx1 95.15408 4.105041 186.2031 45.359624 5.503337
2 xxx2 237.96751 75.001956 400.9331 5.345635 2.418361
pval padj resVarA resVarB
1 7.385318e-139 1.477064e-138 0.07301066 138.9710
2 1.181461e-84 1.181461e-84 5.34503239 266.2948
counts(cds[c("xxx1","xxx2"),])
Control1 Control2 Control3 Contro4 Control5 Control6 Case1
xxx1 2 11 1 4 7 1 7
xxx2 291 51 27 21 109 28 22
Case2 Case3 Case4 Case5 Case6
xxx1 20 1 23 5 990
xxx2 19 15 2298 36 28
This seems surprising, that the p-values should be so small, given the values, barring one outlier, are very similar between control. I suspect the case with a very high value is more likely an outlier, perhaps due to alignment problems or other reasons.
Is there a more robust method to calculate de for RNA seq data?
Many thanks.
Apologies for the bad spacing of the results tables
Just one thing: I'd rather use "padj" instead of "pval"
Absolutely, but even with using padj and running rbinomTest on all genes, the padj goes to e-134
ask Simon Anders, the developper of the package ;)