How to compute RPKM in R? I have the following data of fragment counts for each gene in 16 samples:

&gt; str(expression)
'data.frame':   42412 obs. of  16 variables:
 $ sample1 : int  4555 49 122 351 53 27 1 0 0 2513 ...
 $ sample2 : int  2991 51 55 94 49 10 55 0 0 978 ...
 $ sample3 : int  3762 28 136 321 94 12 15 0 0 2181 ...
 $ sample4 : int  4845 43 193 361 81 48 9 0 0 2883 ...
 $ sample5 : int  2920 24 104 151 50 20 32 0 0 1743 ...
 $ sample6 : int  4157 11 135 324 58 26 4 0 0 2364 ...
 $ sample7 : int  3000 19 155 242 57 12 18 2 0 1946 ...
 $ sample8 : int  5644 30 227 504 91 37 11 0 0 2988 ...
 $ sample9 : int  2808 65 247 93 272 38 1108 1 0 1430 ...
 $ sample10: int  2458 37 163 64 150 29 729 2 1 1049 ...
 $ sample11: int  2064 30 123 51 142 23 637 0 0 1169 ...
 $ sample12: int  1945 63 209 40 171 41 688 3 2 749 ...
 $ sample13: int  2015 57 432 82 104 47 948 4 0 1171 ...
 $ sample14: int  2550 54 177 59 201 36 730 0 0 1474 ...
 $ sample15: int  2425 90 279 73 358 34 1052 3 3 1027 ...
 $ sample16: int  2343 56 365 67 161 43 877 3 1 1333 ...


How do I compute RPKM values from these?
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))
}