Easy way to fetch exons from with Exon stable ID / annotation question

Question

[zeelo]

Hi BioStars,

I have a quick question. I was wondering about a straight-forward command-line way to fetch specific exon sequence data from published genomes, and not the entire gene. E.g., from the Honeybee genome, I want to fetch GB53581-RA-E5, GB52073-RA-E3, GB52625-RA-E6. I found a way for the entire genes but I only want the exons. Is there also a command-line way to identify (fetch?) the orthologous exons/genes from other organisms? Lastly, is there a way to programmatically associate the IDs from the honeybee annotation (say GB52073) with the NCBI Gene ID? (=410059).

Any help is greatly appreciated!

Zeelo

[goodez]

The first thing I would try is to find the gtf file and filter it based on the 3rd column for "exon" features:

zcat organism.gtf.gz | awk '$3=="exon"'

With this list of exon features you can grep for certain genes, then additionally grep for different exon numbers. The 3rd and 4th columns give coordinates for the exon which you can use to fetch the sequences.

One way to find orthologous genes would be to do a command line blast against the other organism.

[Benn]

Do you consider R biomaRt a straight-forward command-line way?

[zeelo]

Hi b.nota,

Thanks for the hint. I explored the package and tried to figure out how it works. However, I couldn't really figure out how to add the ensembl metazoa mart to the package. I guess I must miss something obvious.

Z

Answer 1

[Benn]

A small example with biomaRt in R.

library("biomaRt")

listMarts(host="metazoa.ensembl.org")
             biomart                       version
1       metazoa_mart      Ensembl Metazoa Genes 39
2 metazoa_variations Ensembl Metazoa Variations 39

ensembl <- useMart("metazoa_mart", dataset="amellifera_eg_gene", host="metazoa.ensembl.org")

filters <- listFilters(ensembl)
# type filters to see all

exon_seq <- getBM(attributes=c("ensembl_exon_id", "chromosome_name", "gene_exon"), filters = "ensembl_exon_id", values="GB53581-RA-E5", mart=ensembl, bmHeader=TRUE)
exon_seq
  Chromosome/scaffold name
1                        6
                                                                                                                                                                                                                                                                                                                                                              Exon sequences
1 GATTAATGTAAAATGTTACAAACAAATACATTCATTAACAAGTAAAACATCATTATAATTTGAATTTAATAATAAATACAAAATAATTTTATGAATTTTATGTATTTTACTAATTTTAATATTTTACAAATCCAGCCATTAACATGCTGAAGAGTACAAGTCTTTAATTGTGCTTTCTTTTCACCTAATGAACAACTTCCACCTATGCAATGTCTCCATCTTGTTAGTTTCCCTATACCACAAGTAACACTACAACTTGACCAGTTTCCCCATTGATCCCATATTTTTGGTCCTTGACGACTCATTCTATCTATTATTAAATTATTTGATTTCTAAATAAAAATATTATTACTGGAAGGGTT
  Exon stable ID
1  GB53581-RA-E5

You can of course remove the "chromosome_name", it was just to illustrate that you can add more attributes.

[zeelo]

Thanks a lot, this is great!

[Benn]

I am glad it can help you! If you want orthologs or other IDs check the attributes with listAttributes(ensembl), and make a new query with those attributes.