Calculating RPKM value manually

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

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

[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, 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]

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]

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

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

[TiborNagy]

The calculation looks OK.

[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