How to compute RPKM in R?
First off,

Don’t use RPKMs.

They are truly deprecated because they’re confusing once it comes to paired-end reads. If anything, use FPKMs, which are mathematically the same but use a more correct name (do we count paired reads separately? No, we count fragments).

Even better, use TPM (= transcripts per million), or an appropriate cross-library normalisation method. TMP is defined as:

$$
\text{TPM}_\color{orchid}i =
    {\color{dodgerblue}{\frac{x_\color{orchid}i}{{l_\text{eff}}_\color{orchid}i}}}
    \cdot
    \frac{1}{\sum_\color{tomato}j \color{dodgerblue}{\frac{x_\color{tomato}j}{{l_\text{eff}}_\color{tomato}j}}}
    \cdot
    \color{darkcyan}{10^6}
$$

where


$\color{orchid}i$: transcript index,
$x_i$: transcript raw count,
$\color{tomato}j$ iterates over all (known) transcripts,
$\color{dodgerblue}{\frac{x_k}{{l_\text{eff}}_k}}$: rate of fragment coverage per nucleobase ($l_\text{eff}$ being the effective length),
$\color{darkcyan}{10^6}$: scaling factor (= “per millions”).


That said, FPKM an be calculated in R as follows. Note that most of the calculation happens in log transformed number space, to avoid numerical instability:

fpkm = function (counts, effective_lengths) {
    exp(log(counts) - log(effective_lengths) - log(sum(counts)) + log(1E9))
}


Here, the effective length is the transcript length minus the mean fragment length plus 1; that is, all the possible positions of an average fragment inside the transcript, which equals the number of all distinct fragments that can be sampled from a transcript.

This function handles one library at a time. I (and others) argue that this is the way functions should be written. If you want to apply the code to multiple libraries, nothing is easier using ‹dplyr›:

tidy_expression = tidy_expression %&gt;%
    group_by(Sample) %&gt;%
    mutate(FPKM = fpkm(Count, col_data$Lengths))


However, the data in the question isn’t in tidy data format, so we first need to transform it accordingly using ‹tidyr›:

tidy_expression = expression %&gt;%
    gather(Sample, Count)




This equation fails if all your counts are zero; instead of zeros you will get a vector of NaNs. You might want to account for that.



And I mentioned that TPMs are superior, so here’s their function as well:

tpm = function (counts, effective_lengths) {
    rate = log(counts) - log(effective_lengths)
    exp(rate - log(sum(exp(rate))) + log(1E6))
}