How to define the gene length for RPKM calculation

Question

[Bnf83]

Hi guys,

I would like to calculate the RPKM of my RNA seq experiment. To do this, as from the formula, I need to know the gene length. My starting point are the row reads (single end) counts resulting from:

htseq-count -t exon -i gene_id my_file.sam gencode_v19.gtf > my_count.txt

To calculate the gene length I used the .gtf file from http://www.gencodegenes.org/releases/19.html.

The .gtf file contains start and end positions of genes annotated as: exons, transcript, CDS, start codon etc. I think I have to consider only exons for each gene in order to calculate the length of a gene but I'm not sure. My final purpose is to calculate the gene-level expression of my samples. I'm not interested in splicing variants for example.

Is it right to take into account exons to calculate the gen length for my purpose?

Thank you in advance and kind regards

B

[Benn]

If you use featureCounts from subread (http://subread.sourceforge.net/) instead of htseq for your quantification, you will get the gene length automatically. Just a suggestion.

[xizhesun]

The length in featureCounts is wrong. Because it includes the length of intron!

Answer 1

[lhdcsu]

Different types of gene lengths can be easily calculated using GTFtools:http://www.genemine.org/gtftools.php

type 1: mean of lengths of isoforms.

type 2: median of lengths of isoforms.

type 3: max of lengths of isoforms.

type 4: length of merged exons of isoforms.

Answer 2

[Nathaniel]

If you decide to ignore the existence of different transcripts and work only on the gene metalevel, then do as you suggest: overlap all exons for each gene and sum the length of each exon to obtain an estimate of the gene length

However, as Michael mentions, RPKMs are not the best choice to do differential analysis in gene expression. Take a look at the following video: www.youtube.com/watch?v=TTUrtCY2k-w

P.D. I assume your RNAseq protocol contained the polyA extraction step, to obtain only mature transcripts

Answer 3

[WouterDeCoster]

Since your final purpose is to calculate the gene-level expression I would suggest to go with these counts straight to DESeq2 (or edgeR, but I'm more familiar with DESeq2). You can export normalized counts using far more sophisticated models than RPKM.

Answer 4

[Michael]

There have been several post showing that RPKM can be detrimental instead of doing any good for gene based analysis. For gene DE you wouldn't need it anyway, because gene length doesn't change. Most packages that compute some aggregate (TPM, CPM, FPKM) will compute the gene length for you based on the annotation.

If you want to calculate it yourself approximately, you have already committed to 'exon length' by the use of the parameter. Any read not overlapping exons has no chance of being counted. Some exons might overlap, so it is not good enough to sum all exon lengths, because positions might count multiple times. I would try something like non-overlapping exon length, that means: compute the Union of all exon intervals and sum the length.

[Vasu]

Dear Michael,

I have to calculate rpkm from counts data which I got from STAR. I'm using edgeR for this and it needs gene_length. I'm doing this rpkm calculation not for DE analysis. Could you please tell me how can I calculate rpkm from counts.

Do you think this is the right way to calculate gene_length?

library(GenomicFeatures)
hg19.ens <- makeTxDbFromUCSC(genome="hg19", tablename="ensGene")
exonic <- exonsBy(hg19.ens, by="gene")
red.exonic <- reduce(exonic)
exon.lengths <- sum(width(red.exonic))

[inesdesantiago]

Almost. The last part needs a bit of work in order to work:

library(GenomicFeatures)
hg19.ens <- makeTxDbFromUCSC(genome="hg19", tablename="ensGene")
exonic <- exonsBy(hg19.ens, by="gene")
red.exonic <- reduce(exonic)
exon.lengths <- vapply(width(red.exonic), sum, numeric(1))

Alternatively, if you have an external gtf file , you could do:

library(GenomicFeatures)
txdb <- makeTxDbFromGFF(gtfFile, format="gtf")
exonic <- exonsBy(txdb, by="gene")
red.exonic <- reduce(exonic)
exon.lengths <- vapply(width(red.exonic), sum, numeric(1))