SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Reply
 
Thread Tools
Old 12-05-2011, 06:11 AM   #1
RockChalkJayhawk
Senior Member
 
Location: Rochester, MN

Join Date: Mar 2009
Posts: 191
Default dexseq_prepare_annotation.py being really picky

I'm trying to create a human GTF file using the hg19.gtf file from the TopHat website.The problem is that DEXseq's script to prepare the GTF for DEXSeq is complaining about a couple od different types of errors. These types are the most common:
Code:
Traceback (most recent call last):
  File "dexseq_prepare_annotation.py", line 93, in <module>
    raise ValueError, "Same name found on two strands: %s, %s" % ( str(l[i]), str(l[i+1]) )
ValueError: Same name found on two strands: <GenomicFeature: exonic_part 'FAM95B1' at chr9: 42473935 -> 42474238 (strand '+')>, <GenomicFeature: exonic_part 'FAM95B1' at chr9: 43027964 -> 43027663 (strand '-')>
Code:
Traceback (most recent call last):
  File "dexseq_prepare_annotation.py", line 91, in <module>
    raise ValueError, "Same name found on two chromosomes: %s, %s" % ( str(l[i]), str(l[i+1]) )
ValueError: Same name found on two chromosomes: <GenomicFeature: exonic_part 'YTHDC1' at chr4_ctg9_hap1: 45745 -> 45366 (strand '-')>, <GenomicFeature: exonic_part 'YTHDC1' at chr4: 69180040 -> 69176105 (strand '-')>
Code:
Traceback (most recent call last):
  File "dexseq_prepare_annotation.py", line 89, in <module>
    assert l[i].iv.end <= l[i+1].iv.start, str(l[i+1]) + " starts too early"
AssertionError: <GenomicFeature: exonic_part 'LOC399939' at chr11: 89645253 -> 89644640 (strand '-')> starts too early
What I have been doing so far is to grep out the genes which are causing these issues, but I'm wondering how many more I will have to do by hand to get this file to work. Does anyone already have one?

Here is what I have so far:
Code:
egrep -v "PRAMEF5|LOC728855|LOC646743|YTHDC1|FAM95B1|LOC399939|PRAMEF22|RSPH10B|ANKRD20|chrUn|FLJ20518|RIMBP|LOC100093|GSTT|UGT2A3|SPAG|PMS2L|HIST|LOC440|chr6_|LOC399940|DUX|TRIM|random|chrX|chrY|TMPRSS11E|SNAR|REXO1L|PPP2|LOC727849|AGSK1|FAM25|GOLGA8|SPDYE|RGPD|MIR4283|CDY1|MIR3675|VCY1|GTF2IP|PRY2|FAM41|PPIAL4|SHOX|H2AFB|LOC100288570|LOC440895|TP53TG|DEFB10|LOC100287834|CSAG3|CSF2R|CBWD|LOC728875|GOLGA2P|CTAGE4|NCRNA00230|TISP|LOC642826|RBMY1|OR2A|UGT2B10|XGPY|MIR4650|LOC100133920|MIR3180|PNMA|LOC150527|MIR3179|TTTY|TBC1D|ZNF84|EIF3C|IL3RA|OR4F3|IL9R|LOC100132287|MIR1256|FAM7A2|RNF5P1|CDY2B|MIR1184-1|AGAP9|SSX|CXorf51|LOC100506123|FAM41AY1|RBMY1J|MAGEA2|MIR1244|HSFX|DEFB104B|HIST2H3C|FAM7A|FAM75A|MCART6|chr17_|LIMS3|SPANX|OR4F29|PPIAL4A|ASMTL" hg19.refFlat.gtf >whole_transcriptome.hg19.gtf
RockChalkJayhawk is offline   Reply With Quote
Old 12-05-2011, 06:55 AM   #2
RockChalkJayhawk
Senior Member
 
Location: Rochester, MN

Join Date: Mar 2009
Posts: 191
Default In the mean time

I have a work around that just won't get all the transcripts, but at least the preparation script will work.

Since most of the issues arise because two genes share the same coordinates in the annotation, I first create a list of all gene that overlap on the same strand using the refFlat file:
Code:
cat refFlat.txt |awk '{print $3"\t"$5"\t"$6"\t"$1"\t0\t"$4}'|mergeBed -nms -s -i stdin|grep ";"|cut -f4|perl -ne 's/;/\t/g;split;if (@_[0] ne @_[1]){print}'|perl -pne 's/\t/\n/g'|sort -u>no_de_list.gtf
then, I remove those from the annotation list along with the ones I've already done:
Code:
INPUT_GTF=genes.gtf
fgrep -v -f no_de_list.gtf $INPUT_GTF |egrep -v "FAM138|FOXD4L|OR4F16|UGT2B17|CFC1|UGT2B15|LOC80154|MIR548H3|ANXA8|PLGLB|AMY1|SNORA59|OR7E37P|PRAMEF|LOC728855|LOC646743|YTHDC1|FAM95B1|LOC399939|PRAMEF22|RSPH10B|ANKRD20|chrUn|FLJ20518|RIMBP|LOC100093|GSTT|UGT2A3|SPAG|PMS2L|HIST|LOC440|chr6_|LOC399940|DUX|TRIM|random|chrX|chrY|TMPRSS11E|SNAR|REXO1L|PPP2|LOC727849|AGSK1|FAM25|GOLGA8|SPDYE|RGPD|MIR4283|CDY1|MIR3675|VCY1|GTF2IP|PRY2|FAM41|PPIAL4|SHOX|H2AFB|LOC100288570|LOC440895|TP53TG|DEFB10|LOC100287834|CSAG3|CSF2R|CBWD|LOC728875|GOLGA2P|CTAGE4|NCRNA00230|TISP|LOC642826|RBMY1|OR2A|UGT2B10|XGPY|MIR4650|LOC100133920|MIR3180|PNMA|LOC150527|MIR3179|TTTY|TBC1D|ZNF84|EIF3C|IL3RA|OR4F3|IL9R|LOC100132287|MIR1256|FAM7A2|RNF5P1|CDY2B|MIR1184-1|AGAP9|SSX|CXorf51|LOC100506123|FAM41AY1|RBMY1J|MAGEA2|MIR1244|HSFX|DEFB104B|HIST2H3C|FAM7A|FAM75A|MCART6|chr17_|LIMS3|SPANX|OR4F29|PPIAL4A|ASMTL">whole_transcriptome.hg19.gtf
RockChalkJayhawk is offline   Reply With Quote
Old 12-05-2011, 07:48 AM   #3
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

As a loyal employee of EMBL, I usually download my annotation files from Ensembl rather than from UCSC, and this is why I have tested dexseq_prepare_annotation.gtf only with Ensembl GTF files so far. Hence, if you use the file from here you can get started right away.

Eventually, of course, we will need to put some effort into ensuring that DEXSeq works nicely with files from other sources as well. However, I think that HTSeq's complaints are not completely unsubstantiated. How can a gene have exons on both strands? This is a bit strange, isn't it? Probably, we should add an option to cause the script to simply silently skip such genes.
Simon Anders is offline   Reply With Quote
Old 05-07-2012, 08:38 AM   #4
FuzzyCoder
Member
 
Location: Minnesota

Join Date: Aug 2011
Posts: 13
Default

Simon-

We mapped and aligned our reads using GSNAP and UCSC genome files (actually, Illumina iGenome), and processed the SAM files to make them compatible with DEXSeq's counting script (e.g., stripped 'chr', replaced '.' with 'N' in sequences, etc.).

Is it safe to use an Esembl GTF -> GFF when everything up to this point has been using UCSC? Are there other changes that would need to be made for the two to play nicely together?
__________________
Best Regards,

Paul Bergmann
FuzzyCoder 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 10:39 PM.


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