Hello all, 

I want to calculate RPKM values of my data and, following previous posts, I use the function rpkm() of edgeR. Initially, I checked how the function works on the hypothetical data of http://blog.nextgenetics.net/?e=51 (Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Wagner GP, Kin K, Lynch VJ. Theory Biosci. 2012) and I got as output the "inconsistent" values presented at the second table of  "Inconsistency with RPKM" paragraph of the above webpage. Could someone please advice if there is actually a problem with the rpkm() function in edgeR? if yes, do I have to recalculate the values manually or is there an updated function?

Thank you

Mike

In edgeR, you should run calcNormFactors() before running rpkm(), for example:

y <- DGEList(counts=counts,genes=data.frame(Length=GeneLength))
y <- calcNormFactors(y)
RPKM <- rpkm(y)

Then rpkm will use the normalized effective library sizes to compute rpkm instead of the raw library sizes. This solves the problem pointed out by Wagner et al. 

We view the edgeR approach as better than either raw rpkm or the fix suggested by the paper that you cite.

If you try it out, note though calcNormFactors() is designed to work on real data sets with many genes. It won't necessarily give good results on a toy hypothetical dataset of just a few genes.

There is no problem with the rpkm function in edgeR. It does exactly what it says on the tin, i.e., it computes the reads per kilobase per million for each gene in each sample. The link you provided suggests an adjustment to the RPKMs to avoid the problem of "inconsistency" between samples, but these adjusted values are not RPKMs anymore. If you want this adjustment, you'll just have to do it yourself:

t(t(rpkms)/colSums(rpkms))

for a matrix of rpkms. Personally, I think that these adjusted RPKMs are more difficult to interpret. Consider the example below:

If you compared RPKMs directly between samples A and B, genes 1 and 2 will not be DE (which is the correct state of affairs). However, if you performed the adjustment, you would divide all RPKM values in sample A by 83333333, and those in sample B by 133333333. This would introduce a spurious difference of 60% between A and B for genes 1 and 2, which is not ideal.