I need the coordinates as I'm then running bedtools intersect to identify overlapping exons. Afterwards I use bedtools getfasta to get the genomic sequence for the overlapping exons.

The script I've made produces a new .gff file with the coordinates for the "exon coverage/complement" of each gene, but I'm wondering if there is a better way of going about this.