Calculating RPKM value manually

[ThePresident]

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

[TiborNagy]

The calculation looks OK.

[ThePresident]

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.

[AdrianP]

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

[ThePresident]

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.

[AdrianP]

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.

[ThePresident]

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!

[AdrianP]

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

head file.fastq

what do you get?

[dpryan]

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

[ThePresident]

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

[ThePresident]

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

[Simon Anders]

[deleted] .

[ThePresident]

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

[AdrianP]

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?

[ThePresident]

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.