Differential expression analysis on RPKMs - Contrasts and contrasts of contrasts

[giorgifm]

Dear all,

I'm looking for a tool (better would be a Bioconductor package) able to perform differential expression analysis on a table of RPKMs. I've seen so far only solutions implementing raw read counts, like DEGSeq, edgeR and BaySeq. And yeah, in this particular case I cannot transform RPKMs back to read counts.

Do you have any idea? Furthermore, a great advantage would be the capability of calculating significance of "contrasts of contrasts" analysis, i.e. of tetrafactorial designs.

E.g. having 4 conditions: WT control, WT treated, mutant control, mutant treated.

DE1 = WT treated vs. WT control
DE2 = mutant treated vs. mutant control

DE1DE2 = (mutant trated vs. mutant control) vs. (WT treated vs. WT control)
Something like the great limma does for microarrays. But on RNASeq RPKMs.

Thanks a lot for any hint!

[chadn737]

Why can't you extract the raw read counts?

[peer.b]

It's unfortunately an issue sometimes, since in the beginning RPKM (Mortazavi et al., 2008) were considered a standard and some softwares (e.g. Cufflinks) by default output only FPKM and pairwise condition comparisons. There's an open (more like a closed now) debate on whether one should use RPKM or raw read counts for differential analysis, the latter being the prime choice for statistical reason (see for example the manuals of edgeR and the wiki of this website).
As for your question, given only a table of RPKM, you can try the NOISeq R script (http://bioinfo.cipf.es/noiseq/doku.php?id=tutorial), if you can manage to convince it that your data is already rpkm-normalized. Good luck!

[Simon Anders]

Any valid statistical analysis which only get RPKM values as input will have very little power, because lacking information on whether a certain count value stems from a long gene with many counts or a short gene with few counts, it is bound to assume that all numbers stem from low-count and hence high-variance genes.

I'd be careful with NOISeq, by the way. In the paper on NOISeq, the method is evaluated only on dataset that contain no replicates or only technical replicates, i.e., it is not shown that this method is able to correctly assess biological variability and hence control type-I error rate in a biologically meaningful way.

[dietmar13]

in my hands NOIseq on RAW read counts failed to get significant DE genes called (i.e. only VERY few called), similar to tophat/cuffdiff (only 2 DE genes called) if comparing 12 normal vs 12 colon cancer samples. DESeq, edgeR, BaySeq, SAMseq worked much better...

with edgeR and limma (design matrix) you can analyse your design (wt-mock, wt-treat, mut-mock, mut-treat), but you need raw reads...

[giorgifm]

Thanks everyone for your super-informative answers. I guess I will have to convince my collaborators to move towards raw read counts then.

[Richard Fox]

Check out a recent comparison of methods:

Nookaew et. al. (2012), Nucleic Acids Research, Vol. 40, 10084-10097.

NOISeq does not lack power and compares well with other methods as well as microarrays.

“Although all the methods were in agreement in identifying significant GO terms related to growth (a consequence of the increased specific growth rate during batch cultivations), GO terms known to be relevant during fully respiratory growth are not all in agreement with the different methods. Specifically, edgeR showed some inconsistencies in capturing GO terms associated with fatty acid beta-oxidation terms (as well as DESeq), fatty acid metabolic process and TCA cycle, whereas baySeq weakly identify increased expression of ATP-coupled proton transport and ion transport. Interestingly, the results derived from NOISeq seem to give stronger signals that explain the known differences between batch and chemostat growth better than the results derived from the other methods.”

The authors speculate that edgeR suffers from increased Type I error.