SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Images of 314 chips part 1 jgag0627 Ion Torrent 13 04-06-2016 05:55 PM
PacBio Images part 1 jgag0627 Pacific Biosciences 5 12-24-2012 11:39 PM
Fastq sorting and common part? stoker Bioinformatics 0 07-07-2011 01:53 AM
Sequencing the same part: how many possible way? IgorFobia General 3 03-04-2011 09:40 PM
cufflinks exonic reads count repinementer Bioinformatics 0 08-13-2010 08:25 AM

Reply
 
Thread Tools
Old 02-06-2013, 08:25 AM   #1
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default [DEXSeq] prepare_annotation.py: exonic part starts too early!

When I tried the command below,

"python dexseq_prepare_annotation.py Mt.gtf Mt_dexseq.gtf"

an error occurs, saying that "AssertionError: <GenomicFeature: exonic_part 'Medtr7g090810' at Mt:3458:Mt3.5.1Chr7: 28524762 -> 28523495 (strand '-') > starts too early"

So, would you please give me a tip about how I will be able to fix this error?
Thank you.
syintel87 is offline   Reply With Quote
Old 02-06-2013, 08:58 AM   #2
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Could you grep for 'Medtr7g090810' in you GTF file and post the lines containing the term?
Simon Anders is offline   Reply With Quote
Old 02-06-2013, 09:25 AM   #3
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default

Quote:
Originally Posted by Simon Anders View Post
Could you grep for 'Medtr7g090810' in you GTF file and post the lines containing the term?
The attached file includes the part of 'Medtr7g090810'.
Thank you for sparing your precious time.
Attached Images
File Type: jpg Mt_grep.jpg (80.1 KB, 11 views)
syintel87 is offline   Reply With Quote
Old 02-06-2013, 09:27 AM   #4
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

It seems HTSeq got confused because the same gene occurs on both the "+" and the "-_ strand.
Simon Anders is offline   Reply With Quote
Old 02-06-2013, 09:52 AM   #5
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default

Quote:
Originally Posted by Simon Anders View Post
It seems HTSeq got confused because the same gene occurs on both the "+" and the "-_ strand.
1. So, annotation file has to be fixed??

2. When I have just made small changes on the "dexseq_prepare_annotation.py", it worked.

exons = HTSeq.GenomicArrayOfSets( "auto", stranded=True )
for f in HTSeq.GFF_Reader( gtf_file ):
if f.type != "exon":
continue
f.attr['transcript_id'] = f.attr['transcript_id'].replace( ":", "_" )
exons[f.iv] += ( f.attr['transcript_id'], f.attr['transcript_id'] )

But, seeing the original gtf file, since there are both exon and CDS, I am not sure whether this code is okay for my gtf file or not.

3. I have another gtf file. For this, I also made a small change. And it worked.

exons = HTSeq.GenomicArrayOfSets( "auto", stranded=True )
for f in HTSeq.GFF_Reader( gtf_file ):
if f.type != "CDS":
continue
f.attr['transcript_id'] = f.attr['transcript_id'].replace( ":", "_" )
CDS[f.iv] += ( f.attr['transcript_id'], f.attr['transcript_id'] )

But, I am not sure whether this code is okay or not.

4. Do you think the codes that I have modified would be okay?
The attached file is about original gtf and dexseq_prepare_annotation.py and output gtf.

Thank you very much!
syintel87 is offline   Reply With Quote
Old 02-06-2013, 09:53 AM   #6
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default

Quote:
Originally Posted by Simon Anders View Post
It seems HTSeq got confused because the same gene occurs on both the "+" and the "-_ strand.
This is the attachment.
Attached Files
File Type: txt dexseq_edit.txt (5.9 KB, 13 views)
syintel87 is offline   Reply With Quote
Old 02-06-2013, 10:48 AM   #7
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

No, yopu cannot change from gene ID to transcript ID, because there may be many genes with several overlapping transcripts, and they won't be handled correctly anymore.

You really should fix your GTF file: Wherever the same gene ID is used for features on different strands, add something to the gene ID. If this is complicated, just add a "+" or "-" to all gene IDs.

BTW, dexseq_prepare only looks at "exon" lines and ignored "CDS" lines
Simon Anders is offline   Reply With Quote
Old 02-06-2013, 11:08 AM   #8
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default

Quote:
Originally Posted by Simon Anders View Post
No, yopu cannot change from gene ID to transcript ID, because there may be many genes with several overlapping transcripts, and they won't be handled correctly anymore.

You really should fix your GTF file: Wherever the same gene ID is used for features on different strands, add something to the gene ID. If this is complicated, just add a "+" or "-" to all gene IDs.

BTW, dexseq_prepare only looks at "exon" lines and ignored "CDS" lines
1.
So, you mean I will have to replace each +/- into + ?
(Alternatively, replace each +/- into -).

2.
In my Mhapla.gtf file, it has only exon. So do I need to fix my gtf file by replaceing CDS with exon?
syintel87 is offline   Reply With Quote
Old 02-06-2013, 11:19 AM   #9
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

1. No, change the gene ID from, say, "Medtr7g090810" to "Medtr7g090810+" and "Medtr7g090810-", depending on strand. This is assuming that you know a scripting language. I wouldn't want to do that manually.

Where did you get this strange GTF file from, anyway? Having the same gene name on both strands is a bug.

2. No, why?
Simon Anders is offline   Reply With Quote
Old 02-06-2013, 01:40 PM   #10
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default

Quote:
Originally Posted by Simon Anders View Post
1. No, change the gene ID from, say, "Medtr7g090810" to "Medtr7g090810+" and "Medtr7g090810-", depending on strand. This is assuming that you know a scripting language. I wouldn't want to do that manually.

Where did you get this strange GTF file from, anyway? Having the same gene name on both strands is a bug.

2. No, why?
1. Thank you. I'll try that. I obtained these gtf files from a member of my project group.

2. If Mhapla.gtf only has CDS but dexseq_prepare.py only looks at "exon" lines and ignored "CDS" lines, the output might have no line at all, I guess.
syintel87 is offline   Reply With Quote
Old 02-06-2013, 01:49 PM   #11
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

2. Sure, if it's this way round, you need to change the CDS lines to exon lines.
Simon Anders is offline   Reply With Quote
Old 11-06-2019, 11:38 PM   #12
arkanion
Junior Member
 
Location: Singapore

Join Date: Jul 2016
Posts: 9
Default

Quote:
Originally Posted by Simon Anders View Post
1. No, change the gene ID from, say, "Medtr7g090810" to "Medtr7g090810+" and "Medtr7g090810-", depending on strand. This is assuming that you know a scripting language. I wouldn't want to do that manually.

Where did you get this strange GTF file from, anyway? Having the same gene name on both strands is a bug.

2. No, why?

I have the same problem. If you search UCSC Genome Browser for the gene for ex; HIST2H3C, you will see 2 genes appearing, one on + and on one - strand. dexseq cannot deal with this, but this situation actually happens in reality for those who use gtf file downloaded from UCSC track.
arkanion 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 06:08 PM.


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