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.
3 answers
+1: I've found DESeq to be a pretty robust tool for analysis, and I would always use this step
You could also try the GLM function (the 4th section in the manual, linked below), but I would mostly use that for multivariate analysis (although it works with one variable).
http://bioconductor.org/packages/release/bioc/vignettes/DESeq/inst/doc/DESeq.pdf
If I'm understanding the data correctly, gene xxx1 should not be differentially expressed.
What is happening, I think, is that DESeq is looking at the mean which is artificially increased by the outlier. This is OK but it is also underestimating the variance because it is using an information sharing approach to calculate the variance which makes all the genes have a more "similar" variance. This boosts power when you have 2 or 3 reps, but decreases power when you have more replicates.
You can normalize your data for the number of reads and just use a t-test if you have six reps. The measured variance is the best estimate of variance so it will give you more statistical power and it also does not introduce systematic bias.
A t-test performed better in simulation with 6 reps than DESeq when we looked at it:
http://bioinformatics.oxfordjournals.org/content/29/5/656
in the supplement (as you would expect).
If these numbers are normalized, the p value from a t-test of 0.32 is more correct.
•
written
by
Michele Busby
Log in to answer this question.
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 ;)