2
$\begingroup$

I'm trying to get genome coordinate information programmatically for multiple genes.

As an example, I'd like to get the annotation (list of exons and their chromosomes, starts, and ends) for this UniProt entry: https://www.uniprot.org/uniprotkb/B6TU39/entry.

There is a corresponding GenBank entry, which gives the CDS (https://www.ncbi.nlm.nih.gov/nuccore/EU968504) but I can't see anywhere how to extract the genome coordinates. How can I do this?

Thanks!

$\endgroup$
4
  • $\begingroup$ You have asked for a 'programmatic' solution; what programming language? $\endgroup$
    – gringer
    Commented Oct 15, 2022 at 8:47
  • $\begingroup$ I'll be using Python, but anything using the UniProt/GenBank/etc. REST APIs would work. $\endgroup$ Commented Oct 16, 2022 at 9:21
  • $\begingroup$ In the first instance it'd be useful just to know how to do it using the website GUIs, as that's likely to give enough clues $\endgroup$ Commented Oct 16, 2022 at 9:22
  • $\begingroup$ Talking of web GUIs, you can use Ensembl plants. Put the uniprot accession in the search box and you will find the gene page. $\endgroup$
    – user172818
    Commented Oct 18, 2022 at 17:15

2 Answers 2

2
$\begingroup$

You don't want the EU968504, that's the mRNA not the genome sequence. You want https://www.ncbi.nlm.nih.gov/nuccore/CM007647, also linked to from the UniProt page, just under the link to EU968504:

Screenshot of the UniProt page with the link

You also need the gene name, ZEAMMB73_Zm00001d029604 in this case, available further up on the page in the "Names & Taxonomy" section:

Screenshot of UniProt page showing the relevant section

Next, we can download the GenBank file for the relevant chromosomal contig using NCBI's entrez API:

curl "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=CM007647&rettype=gbwithparts&retmode=text" > CM007647

This will take a few minutes. Once that is done, you can parse the file and extract the start and end coordinates using something like this:

#!/usr/bin/env python3
from Bio import SeqIO
import sys

gb_record=SeqIO.read(open(sys.argv[1],"r"), "genbank")

for feature in gb_record.features:
    if feature.type == "gene" and sys.argv[2] in feature.qualifiers['locus_tag']:
        if feature.location.strand == 1:
            strand = "+"
        elif feature.location.strand == -1:
            strand = "-"
        else:
            print("Unknown strand: ", feature.location.strand)
            sys.exit(1)
            
        print("%s\t%s\t%s\t%s\t%s\t%s\t%s\t.\t." % (gb_record.id, "GenBank",
                                                    feature.type, feature.location.start,
                                                    feature.location.end, ".", strand))

You pass the GenBank file and gene name as arguments, and this is the result:

$ foo.py CM007647  ZEAMMB73_Zm00001d029604
CM007647.1  GenBank gene    78770080    78773320    .   +   .   .
$\endgroup$
2
  • 1
    $\begingroup$ Thanks so much for the detailed answer! $\endgroup$ Commented Oct 16, 2022 at 17:02
  • $\begingroup$ @MichaelDunne CM007647 comes from GRAMENE-4.0, not the latest NAM-5.0 maize genome. On NAM-5.0, this gene is located at 1:77,982,020-77,986,746. $\endgroup$
    – user172818
    Commented Oct 18, 2022 at 17:05
0
$\begingroup$

You could align back the peptide sequence to the genome, respecting splicing, for example using miniprot from Heng Li.

$\endgroup$
2
  • $\begingroup$ As it’s currently written, your answer is unclear. Please edit to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center. $\endgroup$
    – Community Bot
    Commented Oct 17, 2022 at 13:02
  • 1
    $\begingroup$ Hi and welcome to the site! We like answers to be a bit more comprehensive here. Could you expand your answer to explain how to use miniprot to do what the OP asks for? Also, you would need to explain how to go from a protein accession to its fasta sequence to then use that to align to the genome, as well as how to figure out which genome to align to and how to differentiate between homologs. $\endgroup$
    – terdon
    Commented Oct 18, 2022 at 8:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.