How to prepare a right gff file for DEXSeq

[dadada4ever]

Hi all, this seems to be a naive question and any suggestions would be appreciated.
I followed DEXSeq's vignette and it worked well except for a little bit dirt in the step 'Visualization'. I found that even if a gene lies on a negative strand, the exon number showed on x-axis is still counted from left.
As you can see below:

ENSG00000158125, the leftmost is E001 which means exonic part 1.
After comparison between gff converted from Ensembl gtf I realized it seems like when dexseq_prepare_annotation.py handles with gtf it doesn’t look at how the exon number looks like:
In gtf file it’s ENSG00000158125’s last exon (exon 36) but since the coordinate is in an ascending manner so after conversion it becomes first exonic part of this gene in gff file.
In gtf file:

2       protein_coding  exon    31557187        31558874        .       -       .       gene_id "ENSG00000158125"; transcript_id "ENST00000379416"; exon_number "36"; gene_biotype "protein_coding"; gene_name "XDH"; p_id "P84907"; transcript_name "XDH-001"; tss_id "TSS26386";

Converted gff file:

2 genes.gtf exonic_part 31557187 31558874 . - . transcripts "ENST00000379416"; exonic_part_number "001"; gene_id "ENSG00000158125"

Yes it doesn’t affect the main analysis but in Simon’s paper that’s the exact plot I want, see pic below, gene FBgn0004449 lies on negative strand and the exonic part starts from rightmost. Is there any trick to fix my problem here? Since dexseq_prepare_annotation.py doesn’t have any parameters.


Thank you.

[dadada4ever]

Bump. Can anyone share ideas of this problem?

[areyes]

Exon IDs are usually arbitrary, irrespective on the strand of the gene. If you want to change them as in the paper, you could easily change them within R:

data("pasillaExons", package="pasilla")
ecs <- pasillaExons
minusStrandGenes <- names( which( tapply(as.character(fData(ecs)$strand), as.character(fData(ecs)$geneID), unique ) == "-" ) )
geneids <- fData(ecs)$geneID
exonids <- fData(ecs)$exonID

for(x in minusStrandGenes){
   rows <- as.character(geneids) %in% x
   exonids[rows] <- rev(exonids[rows])
}

fData(ecs)$exonID <- exonids
featureNames(ecs) <- paste(geneIDs(ecs), exonIDs(ecs), sep=":")

[dadada4ever]

Thank you Alejandro. I'll try immediately.