SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Remove overlapping feature from GFF ruggedtextile Bioinformatics 1 06-20-2014 05:23 AM
Removing overlapping genes from annotation for RNAseq read count DRAT Bioinformatics 2 04-11-2014 03:53 AM
How do you get the gene length that htseq-count is using? jebe RNA Sequencing 0 03-07-2013 08:37 AM
Gene reads count from RNA-seq sam file by CIGAR Jerry_Zhao Bioinformatics 8 02-01-2013 11:01 AM
non-overlapping features GFF (BEDTools?) seqeve Bioinformatics 2 11-12-2012 10:44 PM

Reply
 
Thread Tools
Old 07-18-2014, 05:47 AM   #1
tompoes
Junior Member
 
Location: Belgium

Join Date: Jul 2014
Posts: 2
Default count non-overlapping gene length from gff file

Hi

I used htseq-count (using default settings) to count all mapped (single-end) reads in a .bam file.
I would like to convert these raw-read count data to RPKM values. Therefore I need to have the non-overlapping transcript (sum of exons) length from all genes in a gff3 file.

You can find a simple example below:

##gff-version
scaffold_1 test gene 47300 48200 . + . ID=A90;Name=A90;
scaffold_1 test mRNA 47300 48200 . + . ID=A90.1;Parent=A90;Name=A90.1;
scaffold_1 test exon 47300 48200 . + . ID=A90.1.1;Parent=A90.1;
scaffold_1 test CDS 47400 48000 . + 0 ID=CDS:A90.1.1;Parent=A90.1;Name=A90.1;
scaffold_1 test gene 48124 49903 . - . ID=A100;Name=A100;
scaffold_1 test mRNA 48124 49903 . - . ID=A100.1;Parent=A100;Name=A100.1;
scaffold_1 test exon 48124 49903 . - . ID=A100.1.1;Parent=A100.1;
scaffold_1 test three_prime_UTR 48124 48464 . - . ID=three_prime_UTR:A100.1.1;Parent=A100.1;Name=A100.1;
scaffold_1 test CDS 48465 49841 . - 0 ID=CDS:A100.1.1;Parent=A100.1;Name=A100.1;
scaffold_1 test five_prime_UTR 49842 49903 . - . ID=five_prime_UTR:A100.1.1;Parent=A100.1;Name=A100.1;

So non overlapping transcript length of A90 is 825 bp (instead of 901) and A100 is 1704 bp (instead of 1780). Are there any tools available that can calculate such unique transcript length or should a script be written?

Thank you for helping!

Last edited by tompoes; 07-18-2014 at 06:02 AM.
tompoes is offline   Reply With Quote
Old 07-18-2014, 07:35 AM   #2
Wallysb01
Senior Member
 
Location: San Francisco, CA

Join Date: Feb 2011
Posts: 286
Default

If all you’re going for is an accurate FPKM, why not just use cufflinks, which assigns FPKMs based on a most likely explanation of the reads method rather than some non-overlapping business gene length combined with htseq counts?

I love htseq/DESeq2 for DE analysis, but for just getting FPKMs cufflinks works great.
Wallysb01 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 03:35 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO