This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Tutorial: Extract Total Non-Overlapping Exon Length Per Gene With Bioconductor

Hi all,

It took me a while to figure this out so I thought it might be useful to a few other people.

When you have used htseq-count on each of your RNA-seq'ed samples and have combined all of your samples in a read count matrix, you might want to turn the absolute read counts into FPKM at some time. Therefore you need the total amount (maybe TMM/quartile-normalized) of mapped reads for each sample (easy) and the total exonic length for each gene (tricky).

The problem is that a gene can have multiple transcripts where the exons of transcript 1 might be partially overlapping with the exons of transcript2. So simply adding up the length of all exons annotated for a specific gene overestimates the gene size while taking the length of just 1 of the transcripts underestimates. Using the GenomicFeatures package from Bioconductor you can easily handle this issue.

# First, import the GTF-file that you have also used as input for htseq-count
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("yourFile.gtf",format="gtf")
# then collect the exons per gene id
exons.list.per.gene <- exonsBy(txdb,by="gene")
# then for each gene, reduce all the exons to a set of non overlapping exons, calculate their lengths (widths) and sum then
exonic.gene.sizes <- sum(width(reduce(exons.list.per.gene)))
fpkm rna-seq bioconductor
• 43,947 views • link
updated by Ram Ram • written by Irsan Irsan

Can different genes that share the same exons cause similar problem?

• 0 views • link
written by PoGibas PoGibas

Because you reduce overlapping exons for each individual gene seperately, they will contribute to both gene sizes like the way you want

• 0 views • link
written by Irsan Irsan

Thanks a lot, save my time writing a perl script. :)

• 0 views • link
written by yyaobo • yyaobo

Nice to hear that it's useful to others

• 0 views • link
written by Irsan Irsan

How about total non-overlapping introns per gene lengths? Introns are grouped by transcripts.How could I group them by genes? Thanks!

• 0 views • link
written by liux.bio liux.bio

Oef, I have to look that up in the GenomicRanges manual. Their is an introns() accessor, that might bring you somewhere

• 0 views • link
written by Irsan Irsan

Hi Irsan, How can you collect exons per gene_name, instead of gene id? I am using a gtf file Thank you

• 0 views • link
written by hoangtv • hoangtv

Do you mean gene symbol in stead of gene name? Like TP53/RB1/E2F3?

• 0 views • link
written by Irsan Irsan

Do you mean gene symbol in stead of gene name? Like TP53/RB1/E2F3?

• 0 views • link
written by Irsan Irsan

Yes, i am using annotation gtf file from Ensembl. I want to get gene symbol, instead of Ensembl gene ID. Do you know how? Thanks

• 0 views • link
written by hoangtv • hoangtv

Many ensembl gene ids do not have an official HGNC symbol (that is what you are looking for). I think the gene symbols are not even in the gtf-file. So you have to work with Ensembl gene IDs and at the end annotate them with gene symbol with for example bioconductor package org.Hs.eg.db

• 0 views • link
written by Irsan Irsan

you're welcome. maybe you should make your post a comment instead of an answer

• 0 views • link
written by Irsan Irsan

Thanks so much, this saves me a lot of time!

• 0 views • link
written by mararie • mararie

Sorry I have wrote I wrong issue and I do not know how to delete it

• 0 views • link
updated by Ram Ram • written by fusion.slope fusion.slope

I used gtf file for human from here but the makeTranscriptDbFromGFF gave me the following error,

error in evaluating the argument 'x' in selecting a method for function 'unique': Error in subs[[i]] : subscript out of bounds

is there any suggestions? probably a point to another gtf file?

• 0 views • link
updated by Ram Ram • written by Quak Quak
Nice it works for you right now. I still wonder what the problem is with the ensembl 75 gtf file. Have you tried removing the first few comment lines (starting with #) in the gtf?
• 0 views • link
written by Irsan Irsan

I've tried this and it might not be the best idea. When creating a GTF-file, the UCSC table browser fills in the gene_ID field with the transcript_ID, no matter whether the UCSC, RefSeq or Ensembl track is used, which defeats the purpose of your script. I have no idea why UCSC would do this, seems like a very dangerous oversight on their part.

Thanks for the script though, it's saved me a lot of time.

• 2 views • link
written by Matthieu Miossec Matthieu Miossec

The last part:

exonic.gene.sizes <- lapply(exons.list.per.gene,function(x){sum(width(reduce(x)))})

Seems to be taking extraordinarily long to run. Any faster method?

• 0 views • link
written by essell22 • essell22
it takes 30 minutes for ensGene table with one 2.67 Ghz cpu
• 0 views • link
written by Irsan Irsan

sum(width(reduce(exons.list.per.gene))) works directly and should be faster

• 0 views • link
written by pgp.martin • pgp.martin

Hi Irsan,

Thanks for a very helpful chunk of code! It solved two of my problems. First I got the length for calculating FPKM which I need later in my analysis, secondly it gave an idea how to quickly count the number of exons per gene which I also needed.

Just a quick comment. Instead of using lapply I simply run:

exonic.gene.sizes <- as.data.frame(sum(width(reduce(exons.list.per.gene))))

Which took only about 2 seconds and returned a data.frame which I could easly merge with the rest of my data.

• 0 views • link
written by ggusiak • ggusiak

updating comment error resolved after starting R new session.

Hi all,

I'm getting an error in this line:

exonic.gene.sizes <- as.data.frame(sum(width(reduce(exons.list.per.gene))))

Error in h(simpleError(msg, call)) : error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': error in evaluating the argument 'x' in selecting a method for function 'width': argument ".f" is missing, with no default

Did anyone encounter this ?

Best

• 0 views • link
written by MatStat MatStat

Hi,

Did you by any chance resolve this? I'm very interested in getting this code to work as well.

Best, Asta

• 0 views • link
written by asta.tranholm • asta.tranholm

I think this error comes from R running "reduce" from dplyr (if you have it loaded) instead of from GenomicRanges as we want it to do. You can substitute "reduce" with "GenomicRanges::reduce" to fix it. So the command will be:

exonic.gene.sizes <- as.data.frame(sum(width(GenomicRanges::reduce(exons.list.per.gene))))
• 0 views • link
written by Acadia • Acadia

Thank you so much for sharing the code! I am a beginner, just wondering if you could help with how to input the length information into the dds data. Many thanks!

• 1 views • link
written by shenwei1376 • shenwei1376

Using the code above For zebrafish, what you can do is:

#first get the gtf file from, lets say, ensembl
library(GenomicFeatures)
#the function above makeTranscriptDbFromGFF is now have different name, make sure that you type the correct one
txdb <- makeTxDbFromGFF("Danio_rerio.GRCz11.93.gtf",format="gtf")
head(txdb)
View(txdb)
# then collect the exons per gene id
exons.list.per.gene <- exonsBy(txdb,by="gene")
head(exons.list.per.gene)
# then for each gene, reduce all the exons to a set of non overlapping exons, calculate their lengths (widths) and sum then
exonic.gene.sizes <- lapply(exons.list.per.gene,function(x){sum(width(reduce(x)))})
head(exonic.gene.sizes)
#to see them in a table format, you should unlist them
unlist_geneLength<-unlist(exonic.gene.sizes)
write.table(unlist_geneLength,"geneLength2.txt")
• 1 views • link
written by realnewbie • realnewbie

The code above was working two days ago, but did not work yesterday. Here is the thing, I guess, one of the packages I have downloaded before GenomicFeatures prevented the code above to work properly (although there were no warning, package masking). Where do I know this from? Simple, today the code is working :). Please, be careful about this.

• 1 views • link
written by realnewbie • realnewbie

Hi Irsan

Thank you very much for providing us with the results of your studies.

The study I am analyzing is a single-end sequence.

Can I use your code with these conditions?

• 0 views • link
written by Maryam • Maryam

Hi,

To reduce the time taken to upload the gtf file. Install using bioconductor the TxDb libraries:

# Gene lenght in bp
## (TxDb.Mmusculus.UCSC.mm10.knownGene)  Mus musculus
## (TxDb.Hsapiens.UCSC.hg19.knownGene) Homo sapiens
## (TxDb.Hsapiens.UCSC.hg38.knownGene) Homo sapiens
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
# then collect the exons per gene id
exons.list.per.gene <- exonsBy(TxDb.Hsapiens.UCSC.hg38.knownGene,by="gene")
# then for each gene, reduce all the exons to a set of non overlapping exons, calculate their lengths (widths) and sum then
exonic.gene.sizes <- sum(width(reduce(exons.list.per.gene)))
# add to count data the gene length. Only works with ENTREZID
data$gene_length<- exonic.gene.sizes[data$entrezid]

best

• 0 views • link
updated by Ram Ram • written by YORDI SEBASTIÁN • YORDI SEBASTIÁN

Would you recommend using the gencode comprehensive annotations file or the basic annotations file to do this since they produce different results (see here and here)? Here they use the comprehensive one.

• 0 views • link
written by gernophil • gernophil

Hi! your code looks great! I want to calculate the and I got to this post via github page but I am confused. Is this code useful to calculate the length only or the meanFragmentLength?

Thank you!

• 0 views • link
written by camillab. • camillab.

Looks like we're on the same path

• 0 views • link
written by Ram Ram

yeah but I am sure you know how to do it, whereas I am just go more confused :/

• 0 views • link
written by camillab. • camillab.

Nope, I kinda gave up. I got to the image on what fragment size is (What is the difference between a Read and a Fragment in RNA-seq?) but cannot figure out how to connect that to the picard metrics output. For the moment, I'd say use insert size for fragment length especially since you trimmed adapters using trimmomatic anyway.

• 0 views • link
written by Ram Ram

4 answers

Here's another way to do it with Python instead of R:

• 0 views • link
written by Kamil Kamil

Thanks for this. Can I just triple-check that the code "coding_lengths.py" is returning the total length of all of the exons (I am trying to calculate RPKM, so I think this is perfect to calculate gene lengths for that purpose).

And then just as a side note, is this data available for other species? I'm on the FTP site for "ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2ensembl.gz" now, and I can seeing anything resembling a folder for other species? Similar for the Sanger FTP site, I can see human and mouse; my work involves human mouse and rat, I'm wondering if this script would work on all of these species, and if not, if someone could suggest an alternative method?

• 0 views • link
written by Tom • Tom

I updated the comment in count_bp with a visual example of how exons are counted.

You can find gene annotations in GTF format for many species here:

ftp://ftp.ensembl.org/pub/release-87/gtf/

Note that the version (87 at the time of this post) is continually updated over time.

Regarding converting Ensembl ID to Entrez ID or vice versa, I recommend you use Ensembl Biomart.

• 0 views • link
written by Kamil Kamil

Thanks. I actually tried to run the script using exactly as above and I ran into an error. In the attached picture: picture, in the purple section I have shown the first few files of the .gtf file, and in the orange section the gene2ensembl file, just to show that I've just downloaded these exactly as you said. Then the green section I have run the code.

The code is also here:

python coding_lengths.py -g gencode.v19.annotation.gtf -n gene2ensembl -o output

Wed Jan 18 11:30:46 2017 # Reading the Gencode annotation file: gencode.v19.annotation.gtf
Wed Jan 18 11:34:36 2017 # done in 230 s
coding_lengths.py:81: FutureWarning: sort(columns=....) is deprecated, use sort_values(by=.....)
  exon.sort(['seqname','start','end'], inplace=True)
Wed Jan 18 11:34:39 2017 # Calculating coding region (exonic) length for each gene...
Wed Jan 18 11:38:18 2017 # done in 219 s
Wed Jan 18 11:38:18 2017 # Reading NCBI mapping of Entrez GeneID to Ensembl gene identifier: gene2ensembl
Traceback (most recent call last):
  File "coding_lengths.py", line 163, in <module>
    main(args)
  File "coding_lengths.py", line 99, in main
    'Ensembl_protein_identifier'])
  File "/home/nis/aoife/env/local/lib/python2.7/site-packages/pandas/io/parsers.py", line 498, in parser_f
    return _read(filepath_or_buffer, kwds)
  File "/home/nis/aoife/env/local/lib/python2.7/site-packages/pandas/io/parsers.py", line 275, in _read
    parser = TextFileReader(filepath_or_buffer, **kwds)
  File "/home/nis/aoife/env/local/lib/python2.7/site-packages/pandas/io/parsers.py", line 590, in __init__
    self._make_engine(self.engine)
  File "/home/nis/aoife/env/local/lib/python2.7/site-packages/pandas/io/parsers.py", line 731, in _make_engine
    self._engine = CParserWrapper(self.f, **self.options)
  File "/home/nis/aoife/env/local/lib/python2.7/site-packages/pandas/io/parsers.py", line 1103, in __init__
    self._reader = _parser.TextReader(src, **kwds)
  File "pandas/parser.pyx", line 515, in pandas.parser.TextReader.__cinit__ (pandas/parser.c:4948)
  File "pandas/parser.pyx", line 705, in pandas.parser.TextReader._get_header (pandas/parser.c:7386)
  File "pandas/parser.pyx", line 829, in pandas.parser.TextReader._tokenize_rows (pandas/parser.c:8838)
  File "pandas/parser.pyx", line 1833, in pandas.parser.raise_parser_error (pandas/parser.c:22649)
pandas.parser.CParserError: Error tokenizing data. C error: Calling read(nbytes) on source failed. Try engine='python'.

I know it is not your job to solve my system errors. My only question is; to a more experienced coder, does this look like a problem with the script itself, or does it look like a problem on my end (i.e. with the way packages are installed or something)?

Thanks

• 0 views • link
written by Tom • Tom

My script expects to read a file compressed with gzip, but you provided a decompressed file.

Try running gzip gene2ensembl and then run the script again.

• 0 views • link
written by Kamil Kamil

Thank you, you were so helpful, the script is absolutely fantastic.

• 0 views • link
written by Tom • Tom

Just a small note to everyone, the Ensembl GTF files are not the same as the NCBI Gencode GTF files. For this script, the only difference you seem to need to make is in the coding_lengths.py script, change line 73 to:

idx = (gc.feature == 'exon') & (gc.transcript_biotype == 'protein_coding')

(so you're changing gc.transcript_type to gc.transcript.biotype).

• 0 views • link
written by Tom • Tom

Recommend the GTFtools package at: http://www.genemine.org/gtftools.php.

Using the '-l' options, four types of gene lengths will be calculated. One of them is the non-overlapping exon length (merge all exons of multiple isoforms of the same gene).

• 0 views • link
written by lhdcsu • lhdcsu

Thanks for the scripts. But when I run the 2nd line "exons.list.per.gene", I get "error, RS-DBI driver: (expired SQLiteConnection)". Dose anyone know why?

• 0 views • link
updated by Ram Ram • written by yc1123 • yc1123

Hi! in case anyone ends here looking to calculate non-overlapping exon lengths, I recently develop this: https://github.com/alejandrogzi/noel, an extremely fast non-overlapping exon length calculator written in Rust.

This program outperforms all the script/tools already mentioned in this discussion, check the tool post here: NOEL: An extremely fast Non-Overlapping Exon Length calculator written in Rust.

• 0 views • link
written by alejandrogzi • alejandrogzi

Log in to answer this question.

More posts like this