Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Calculating RPKM value manually

    Hello,

    I'm dealing with bacterial RNA-seq data. I would like to show absolute gene expression by calculating RPKM value for every gene. Initially, I've been using Cufflinks for that, but I don't like the way Cufflinks deals with it, so I decided to manually calculate my iwn RPKM values.

    The equation I want to use is:

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

    C = Number of reads mapped to a gene
    N = Total mapped reads in the experiment
    L = gene length in base-pairs for a gene

    Now, my question is the following: I want to use raw counts (obtained from HTSeq count table) as the number of reads that actually mapped to a gene (C). I think that should be okay. However, for the total mapped reads in the experiment (N) I thought about simply adding all reads from the HTSeq table.

    For exemple:

    If my gene length (L) is : 200pb
    Number of reads mapped (as from HTSeq-count) (C) : 400
    Total mapped reads (sum for all genes from HTSeq-count) (N) : 10^8

    RPKM = (10^9 * 400)/(10^8 * 200) = 20

    Would that be the right way of calculating it? My aim is to do the same thing for each of my sequencing libraries and then simply compare.

    Thanks you guys,

    TP

  • #2
    The calculation looks OK.

    Comment


    • #3
      Yeah, the more I think about it, the more it makes sense. HTSeq-count uses alignment files from bowtie and DESeq (used for DE analysis) uses HTSeq-count tables... Thus, there seems to be no reason why I couldn't use raw reads from HTSeq-count table for RPKM calculations.

      Comment


      • #4
        Is your data paired end? In that case, you need FPKM counts, not RPKM.

        Comment


        • #5
          Originally posted by AdrianP View Post
          Is your data paired end? In that case, you need FPKM counts, not RPKM.
          Nope, single-end... That's why I don't like using Cufflinks who gives FPKM values.

          Comment


          • #6
            Originally posted by ThePresident View Post
            Nope, single-end... That's why I don't like using Cufflinks who gives FPKM values.
            In case of single-end reads, FPKM=RPKM. But check and see if the values match, I am curious.

            Comment


            • #7
              Nope, it does not match for one obvious reason (and I'm again banging my head against the wall):

              - I can't just sum all the reads from HTseq table: the number is way over the total number of reads from the library. I'm trying to see why... damn!

              Comment


              • #8
                Originally posted by ThePresident View Post
                Nope, it does not match for one obvious reason (and I'm again banging my head against the wall):

                - I can't just sum all the reads from HTseq table: the number is way over the total number of reads from the library. I'm trying to see why... damn!
                Use the .fastq file with all of your reads. When you do
                Code:
                head file.fastq
                what do you get?

                Comment


                • #9
                  Are you counting multimappers? Also, cufflinks calulated FPKMs will almost never be the same as those calculated by hand (partly due to using a different gene length and partly due to cufflinks using fractional read counts).

                  Comment


                  • #10
                    @AdrianP : I'm at work now and I don't have access to my linux station now. But, from bowtie alignment I know how much reads I have in my libraries (fastq). For example, in one of my libraries I have 61,402,323 reads, and when counted with HTSeq-count (sum of all reads across all genes) I get 116,898,233 which is almost double.

                    @dpryan: My understanding of HTSeq-count is that it does not count multimappers. I used intersection-nonempty mode. Here is the final output:

                    no_feature 9034352
                    ambiguous 958299
                    too_low_aQual 0
                    not_aligned 2097930
                    alignment_not_unique 0

                    Fastq file contained 61,402,323 reads and 96.58% of that number have been aligned with bowtie. So, the SAM file has to contain around 59,302,363. So, if there are multimapped reads in mu HTSeq-count, how to get rid of them?

                    Comment


                    • #11
                      But on the other side, the total number of mapped reads in a the RPKM formula is quite arbitrary. It will be the same for every gene inside one single library, so as long as keep that in mind, the formula still makes sense.

                      On the other side, I wonder if I could compare RPKM values calculated by this manner in between libraries containing replicates? Ex. If I have three libraries of the same condition (biological replicates) but with different sequencing depth, could I calculate RPKM as explained above and then just average them...?

                      Comment


                      • #12
                        [deleted] .
                        Last edited by Simon Anders; 01-24-2014, 10:36 AM. Reason: post was meant to appear in another thread, but I confused my browser tabs

                        Comment


                        • #13
                          It might seems that I'm flooding this thread, but I just realized that I already had the same problem. I think that the right way of dealing with all this is to extract only uniquely-mapped reads from bowtie-generated BAM files and use them in HTSeq-count, DESeq and RPKM calculations.

                          My only concern is to lose too much reads du to the multimapping...

                          Comment


                          • #14
                            Originally posted by ThePresident View Post
                            It might seems that I'm flooding this thread, but I just realized that I already had the same problem. I think that the right way of dealing with all this is to extract only uniquely-mapped reads from bowtie-generated BAM files and use them in HTSeq-count, DESeq and RPKM calculations.

                            My only concern is to lose too much reads du to the multimapping...
                            If you use 61,402,323 as total reads, what value of RPKM do you get for any given gene compared to FPKM by cufflinks?

                            Comment


                            • #15
                              Originally posted by AdrianP View Post
                              If you use 61,402,323 as total reads, what value of RPKM do you get for any given gene compared to FPKM by cufflinks?
                              It's somewhat near but not exactly the same thing. I did it with just a couple of values but still: 8vs9 , 146vs300, 13vs18.

                              Comment

                              Latest Articles

                              Collapse

                              • seqadmin
                                Essential Discoveries and Tools in Epitranscriptomics
                                by seqadmin


                                The field of epigenetics has traditionally concentrated more on DNA and how changes like methylation and phosphorylation of histones impact gene expression and regulation. However, our increased understanding of RNA modifications and their importance in cellular processes has led to a rise in epitranscriptomics research. “Epitranscriptomics brings together the concepts of epigenetics and gene expression,” explained Adrien Leger, PhD, Principal Research Scientist on Modified Bases...
                                Yesterday, 07:01 AM
                              • seqadmin
                                Current Approaches to Protein Sequencing
                                by seqadmin


                                Proteins are often described as the workhorses of the cell, and identifying their sequences is key to understanding their role in biological processes and disease. Currently, the most common technique used to determine protein sequences is mass spectrometry. While still a valuable tool, mass spectrometry faces several limitations and requires a highly experienced scientist familiar with the equipment to operate it. Additionally, other proteomic methods, like affinity assays, are constrained...
                                04-04-2024, 04:25 PM

                              ad_right_rmr

                              Collapse

                              News

                              Collapse

                              Topics Statistics Last Post
                              Started by seqadmin, 04-11-2024, 12:08 PM
                              0 responses
                              37 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 04-10-2024, 10:19 PM
                              0 responses
                              41 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 04-10-2024, 09:21 AM
                              0 responses
                              35 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 04-04-2024, 09:00 AM
                              0 responses
                              54 views
                              0 likes
                              Last Post seqadmin  
                              Working...
                              X