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 %>%
group_by(Sample) %>%
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 %>%
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))
}
