# How to compute RPKM in R?

I have the following data of fragment counts for each gene in 16 samples:

> 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?

• What did you tried to solve this question? Did you visit Bioconductor?
– llrs
May 17 '17 at 10:46

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))
}

• Can I ask sort of a meta question? I've started seeing '%>%' in R code recently and had never noticed it before. What does that do exactly?
– Greg
May 17 '17 at 18:36
• @Greg It’s a pipe. It has been around for quite a long time (though different libraries used different operator symbols; I myself used to use %|%, similar to the shell pipe) but only recently gained mainstream traction, mainly through the dplyr library. May 17 '17 at 20:40

RPKM is defined as:

As you can see, you need to have gene lengths for every gene.

Let's say geneLength is a vector which have the same number of rows as your data.frame, and every value of the vector corresponds to a gene (row) in expression.

expression.rpkm <- data.frame(sapply(expression, function(column) 10^9 * column / geneLength / sum(column)))


Regarding numerical stability

It is suggested in one of the answers, that all the computations should be performed in a log-transformed scale. In my opinion there is absolutely no reason for doing that. IEEE binary64 stores a number as binary number 1.b_{51}b{50}...b_0 times 2^{e-1023}. The relative precision doesn't depend on the exponent value given that a number is in the range [~10^-308; 10^308].

In case of RPKM we can only get out of the range if total number of reads is around 10^300, which is not realistic at all.

On the bright site there is not much harm in doing computations in the log-scale either. Worst case you loose a bit of precision.

• and then picking the right gene length itself can be a nightmare depending on the organism... So maybe you should also develop what could be the strategies to pick the correct one. May 17 '17 at 10:33
• @Mitra I think this is slightly out of scope. But if you have an experience in doing that, could you add an answer here? Would be really great! May 17 '17 at 10:41
• @Mitra given it a second thought, this deserves a separate question. If there will be one, I'll link it from this answer. May 17 '17 at 10:44
• I don't think it is out of scope: you need the gene length to compute the RPKM. Other answers bypassed this problem by using transcript or exon length which are easily retrieved from the annotation. That is not the case with gene length which you are using in your formula. May 19 '17 at 13:29

If you are planning to do a differential expression analysis, you will probably don't need the RPKM calculation.

RPK= No.of Mapped reads/ length of transcript in kb (transcript length/1000)

RPKM = RPK/total no.of reads in million (total no of reads/ 1000000)

The whole formula together:

RPKM = (10^9 * C)/(N * L) Where,

C = Number of reads mapped to a gene

N = Total mapped reads in the experiment

L = exon length in base-pairs for a gene

If you are looking for a more visual solution (in addition to the other answers), NCI Genomic Data Commons (TCGA repository) offers a nice formula:

Where:

RCg: Number of reads mapped to the gene

RCpc: Number of reads mapped to all protein-coding genes

RCg75: The 75th percentile read count value for genes in the sample

L: Length of the gene in base pairs

As others have pointed out, FPKMs have some problems. GDC also calculates FPKM-UQ values that are upper quartile normalized. Those are recommended for cross-sample comparison and differential expression analysis.

• The denominator for FPKM-UQ is still sample specific, so this will not be an appropriate normalisation for differential expression analysis (it's within sample normalisation, not between). May 17 '17 at 15:46
• GDC seems to imply otherwise: gdc.cancer.gov/about-data/data-harmonization-and-generation/… May 17 '17 at 17:17