Answers.Body: <p>RPKM is defined as:</p>
  
  <blockquote>
    <p>RPKM =   numberOfReads / ( geneLength/1000 * totalNumReads/1,000,000 )</p>
  </blockquote>
  
  <p>As you can see, you need to have gene lengths for every gene.</p>
  
  <p>Let's say <code>geneLength</code> is a vector which have the same number of rows as your <code>data.frame</code>, and every value of the vector corresponds to a gene (row) in <code>expression</code>.</p>
  
  <pre><code>expression.rpkm &lt;- data.frame(sapply(expression, function(column) 10^9 * column / geneLength / sum(column)))
  </code></pre>
  
  <p><strong>Regarding numerical stability</strong></p>
  
  <p>It is suggested in <a href="https://bioinformatics.stackexchange.com/a/69/191">one of the answers</a>, that all the computations should be performed in a log-transformed scale. In my opinion there is absolutely no reason for doing that. IEEE <a href="https://en.wikipedia.org/wiki/Double-precision_floating-point_format" rel="nofollow noreferrer">binary64</a> stores a number as binary number 1.b_{51}b{50}...b_0 times 2^{e-1023}. The <em>relative</em> precision doesn't depend on the exponent value given that a number is in the range [~10^-308; 10^308].</p>
  
  <p>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.</p>
  
  <p>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.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T10:33:45.303Z
Answers.Text: 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.
Answers.UserDisplayName: Mitra
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T10:41:34.350Z
Answers.Text: @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!
Answers.UserDisplayName: Iakov Davydov
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T10:44:03.900Z
Answers.Text: @Mitra given it a second thought, this deserves a separate question. If there will be one, I'll link it from this answer.
Answers.UserDisplayName: Iakov Davydov
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-19T13:29:29.870Z
Answers.Text: 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.
Answers.UserDisplayName: Mitra
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2017-05-17T10:28:32.203Z
Answers.LastActivityDate: 2018-07-26T14:43:27.470Z
Answers.LastEditDate: 2018-07-26T14:43:27.470Z
Answers.UserDisplayName: Iakov Davydov
Answers.Body: <p>If you are planning to do a differential expression analysis, you will probably don't need the RPKM calculation.</p>
  
  <p>RPK= No.of Mapped reads/ length of transcript in kb (transcript length/1000)</p>
  
  <p>RPKM = RPK/total no.of reads in million (total no of reads/ 1000000)</p>
  
  <p><strong>The whole formula together:</strong></p>
  
  <p>RPKM = (10^9 * C)/(N * L)
  Where,</p>
  
  <p>C = Number of reads mapped to a gene</p>
  
  <p>N = Total mapped reads in the experiment</p>
  
  <p>L = exon length in base-pairs for a gene</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T15:33:34.013Z
Answers.LastActivityDate: 2017-05-17T15:33:34.013Z
Answers.UserDisplayName: arup
Answers.Body: <p>If you are looking for a more visual solution (in addition to the other answers), <a href="https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/" rel="nofollow noreferrer">NCI Genomic Data Commons (TCGA repository)</a> offers a nice formula:</p>
  
  <p><a href="https://i.sstatic.net/RY53Q.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/RY53Q.png" alt="enter image description here"></a></p>
  
  <blockquote>
    <p>Where:</p>
    
    <p>RCg: Number of reads mapped to the gene</p>
    
    <p>RCpc: Number of reads mapped to all protein-coding genes</p>
    
    <p>RCg75: The 75th percentile read count value for genes in the sample</p>
    
    <p>L: Length of the gene in base pairs</p>
  </blockquote>
  
  <p>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.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T15:46:34.070Z
Answers.Text: 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).
Answers.UserDisplayName: sjcockell
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T17:17:09.197Z
Answers.Text: GDC seems to imply otherwise: https://gdc.cancer.gov/about-data/data-harmonization-and-generation/genomic-data-harmonization/high-level-data-generation/rna-seq-quantification
Answers.UserDisplayName: burger
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T15:16:10.923Z
Answers.LastActivityDate: 2017-05-17T15:16:10.923Z
Answers.UserDisplayName: burger
Answers.Body: <p>First off,</p>
  
  <p><strong>Don’t use RPKMs</strong>.</p>
  
  <p><a href="https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/" rel="noreferrer">They are truly deprecated</a> because they’re confusing once it comes to paired-end reads. If anything, use <em>FPKM</em>s, which are mathematically the same but use a more correct name (do we count paired reads separately? No, we count <em>fragments</em>).</p>
  
  <p>Even better, <a href="http://rpubs.com/klmr/rnaseq-norm" rel="noreferrer">use TPM (= transcripts per million), or an appropriate cross-library normalisation method</a>. TMP is defined as:</p>
  
  <p>$$
  \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}
  $$</p>
  
  <p>where</p>
  
  <ul>
  <li>$\color{orchid}i$: transcript index,</li>
  <li>$x_i$: transcript raw count,</li>
  <li>$\color{tomato}j$ iterates over all (known) transcripts,</li>
  <li>$\color{dodgerblue}{\frac{x_k}{{l_\text{eff}}_k}}$: rate of fragment coverage per nucleobase ($l_\text{eff}$ being the <a href="https://bioinformatics.stackexchange.com/q/367/29">effective length</a>),</li>
  <li>$\color{darkcyan}{10^6}$: scaling factor (= “per millions”).</li>
  </ul>
  
  <p>That said, FPKM an be calculated in R as follows. Note that most of the calculation happens in log transformed number space, to avoid <a href="https://en.wikipedia.org/wiki/Numerical_stability" rel="noreferrer">numerical instability</a>:</p>
  
  <pre class="lang-r prettyprint-override"><code>fpkm = function (counts, effective_lengths) {
      exp(log(counts) - log(effective_lengths) - log(sum(counts)) + log(1E9))
  }
  </code></pre>
  
  <p>Here, the <em>effective length</em> 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.</p>
  
  <p>This function handles <em>one</em> library at a time. I (<a href="http://varianceexplained.org/r/tidy-genomics-biobroom/" rel="noreferrer">and others</a>) argue that this is the way functions should be written. If you want to apply the code to multiple libraries, nothing is easier using <a href="http://dplyr.tidyverse.org/" rel="noreferrer">‹dplyr›</a>:</p>
  
  <pre><code>tidy_expression = tidy_expression %&gt;%
      group_by(Sample) %&gt;%
      mutate(FPKM = fpkm(Count, col_data$Lengths))
  </code></pre>
  
  <p>However, the data in the question isn’t in tidy data format, so we first need to transform it accordingly using <a href="http://tidyr.tidyverse.org/" rel="noreferrer">‹tidyr›</a>:</p>
  
  <pre><code>tidy_expression = expression %&gt;%
      gather(Sample, Count)
  </code></pre>
  
  <hr>
  
  <p>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.</p>
  
  <hr>
  
  <p>And I mentioned that TPMs are superior, so here’s their function as well:</p>
  
  <pre class="lang-r prettyprint-override"><code>tpm = function (counts, effective_lengths) {
      rate = log(counts) - log(effective_lengths)
      exp(rate - log(sum(exp(rate))) + log(1E6))
  }
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T18:36:16.750Z
Answers.Text: 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?
Answers.UserDisplayName: Greg
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T20:40:12.327Z
Answers.Text: @Greg [It’s a pipe](https://github.com/tidyverse/magrittr). 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.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T10:46:41.653Z
Answers.LastActivityDate: 2017-07-26T16:13:26.200Z
Answers.LastEditDate: 2017-07-26T16:13:26.200Z
Answers.UserDisplayName: Konrad Rudolph
Body: <p>I have the following data of fragment counts for each gene in 16 samples:</p>
  
  <pre><code>&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 ...
  </code></pre>
  
  <p>How do I compute RPKM values from these?</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-17T10:46:03.897Z
Comments.Text: What did you tried to solve this question? Did you visit Bioconductor?
Comments.UserDisplayName: llrs
Tags: <rna-seq><r><transcriptome><fpkm>
Title: How to compute RPKM in R?
ThreadText: # How to compute RPKM in R?
  
  (Asked by: username_0 on 2017-05-17)
  
  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?
  
  **Tags:** rna-seq, r, transcriptome, fpkm
  
  > **Comments:**
  > **username_1**: What did you tried to solve this question? Did you visit Bioconductor?  
  
  ---
  
  ## 4 answers
  
  ### Answer by username_6 on 2017-05-17. Score: 31 | Accepted answer
  
  First off,
  
  **Don’t use RPKMs**.
  
  [They are truly deprecated](https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/) because they’re confusing once it comes to paired-end reads. If anything, use *FPKM*s, 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](http://rpubs.com/klmr/rnaseq-norm). 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](https://bioinformatics.stackexchange.com/q/367/29)),
  * $\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](https://en.wikipedia.org/wiki/Numerical_stability):
  
  ```r
  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](http://varianceexplained.org/r/tidy-genomics-biobroom/)) 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›](http://dplyr.tidyverse.org/):
  
  ```
  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›](http://tidyr.tidyverse.org/):
  
  ```
  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:
  
  ```r
  tpm = function (counts, effective_lengths) {
      rate = log(counts) - log(effective_lengths)
      exp(rate - log(sum(exp(rate))) + log(1E6))
  }
  
  ```
  
  > **Comments:**
  > **username_7**: 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?  
  > **username_6**: @username_7 [It’s a pipe](https://github.com/tidyverse/magrittr). 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.  
  
  ---
  
  ### Answer by username_0 on 2017-05-17. Score: 8
  
  RPKM is defined as:
  
  > RPKM = numberOfReads / ( geneLength/1000 \* totalNumReads/1,000,000 )
  
  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](https://bioinformatics.stackexchange.com/a/69/191), 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](https://en.wikipedia.org/wiki/Double-precision_floating-point_format) 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.
  
  > **Comments:**
  > **username_2**: 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.  
  > **username_0**: @username_2 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!  
  > **username_0**: @username_2 given it a second thought, this deserves a separate question. If there will be one, I'll link it from this answer.  
  > **username_2**: 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.  
  
  ---
  
  ### Answer by username_3 on 2017-05-17. Score: 1
  
  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
  
  
  ---
  
  ### Answer by username_4 on 2017-05-17. Score: 0
  
  If you are looking for a more visual solution (in addition to the other answers), [NCI Genomic Data Commons (TCGA repository)](https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/) offers a nice formula:
  
  [![enter image description here](https://i.sstatic.net/RY53Q.png)](https://i.sstatic.net/RY53Q.png)
  
  > 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.
  
  > **Comments:**
  > **username_5**: 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).  
  > **username_4**: GDC seems to imply otherwise: https://gdc.cancer.gov/about-data/data-harmonization-and-generation/genomic-data-harmonization/high-level-data-generation/rna-seq-quantification  
  
  ---
  
  