11

I have read counts data and I want to convert them into RPKM values. For this conversion I need the gene length.

Does the gene length need to be calculated based on the sum of coding exonic lengths? Or are there any different ways for that?

I know that gene length can be taken from the Gencode GTF v19 file. Could you please tell me how that Gene_length is calculated?

The Data I'm having is RNA-Seq data. I don't have any idea whether I need to include UTR's in this calculation or only exons?

In Github I have seen RPKM calculation from Counts data with the Gene_length from Gencode GTF file. Do you think this is the right way of calculation?

And why RPKM is - Its not for differential analysis. For TNBC subtyping they use microarray data. I would like to give a try with RNA-Seq data. So for this I'm trying out different and the right way.

Konrad Rudolph
  • 4,845
  • 14
  • 45
stack_learner
  • 1,262
  • 14
  • 26

2 Answers2

7

Here you can find some example R code to compute the gene length given a GTF file (it computes GC content too, which you don't need). This uses one of a number of ways of computing gene length, in this case the length of the "union gene model". In this method, the non-duplicated exons for each gene are simply summed up ("non-duplicated" in that no genomic base is double counted). This is a very simple way of getting a gene length.

There are alternative methods that you should be aware of, among which are:

  1. Median transcript length: That is, the exonic lengths in each transcript are summed and the median across transcripts is used. This is probably a little more valid than the code that I linked to.
  2. Per-sample effective gene lengths: the optimal method, though it requires using something like RSEM, which will give you an effective gene length.
  3. Using the length of the "major isoform" in your tissue of interest. This isn't as good as method 2, but is more accurate than all of the others.

At the end of the day, you're just coming up with a scale factor for each gene, so unless you intend to compare values across genes (this is problematic to begin with) then it's questionable if using some of the more correct but also more time-involved methods are really getting you anything.

Edit: Note that if you want to plug these values into some sort of subtyping tool (TNBC in your case), you should first start with some samples for which you know the subtype. Then you can at least see if you're getting reasonable results. After that, do read up on how the method works and see if there's anything about RNAseq that makes it incompatible.

Devon Ryan
  • 19,602
  • 2
  • 29
  • 60
  • I would think that the method used to calculate gene length should be informed by the counting method. If reads were counted across all exons, does it make much sense to use the alternative methods you mention? – heathobrien Sep 26 '17 at 10:12
  • There are data-dependent methods (namely option 2 and maybe 3) and data-independent methods (everything else). The counting method is irrelevant except with things like RSEM which are going to produce effective lengths based on the relative transcript expression observed in each sample. Otherwise, a gene's length is just a constant. – Devon Ryan Sep 26 '17 at 10:52
1

I assume you are mapping against the genome rather the transcriptome, since for the later the length would be trivial.

Assuming the first, I think not only the coding sections should be included but also the UTR, since reads can map against them which is what we ultimately care about.

In general, I found gene annotation files (e.g. gff or gtf) can be inconsistent in terms of naming, so it's good practice to inspect and double check. Below is some R code to import the annotation and calculate isoform lengths:

library(rtracklayer)
# Reading into a GRanges object
anno <- import.gff3("annotation.gff")
# Filtering exons and UTRs
exons <- anno[anno@elementMetadata$type %in% c("exon","five_prime_UTR","three_prime_UTR"),]

Depending on the annotation at hand, the most sensible is probably best to count the length of each isoform which are often contained in the "Parent" column of the annotation file:

# splitting up isoforms as preparation for the next step
tmp <- split(exons,as.character(exons$Parent))
# for each isoform, calculate the sum of all reduced exons
Gene_length <- sum(width(reduce(tmp)))

Note, reduce merges overlapping intervals together, since UTRs can "contain" bits of exons which would be otherwise double counted.

This code can of course be adapted mainly by changing the "Parent", "exon" etc.

Sebastian Müller
  • 600
  • 2
  • 9
  • If you're filtering for exons then you needn't include the UTRs. You're not hurting anything since you reduce() them out anyway, but you could have freed up some memory. You can also do the filtering directly in import.gff3 (see the first link in my answer). – Devon Ryan Oct 16 '17 at 07:07