I have a genomic region of interest and I am trying to retrieve the gene (if its in one) and amino acid sequence (if its in a translated region). Suppose my region of interest is:

chromosome = "chr19"
start = 55014736
stop = 55014753

If I put this into USCS genome browser by hand I see it's in gene GP6 and the translation is LCVSLC. I can use cruzdb and retrieve the correct gene but the translation of cds_sequence does not contain the LCVSLC motif. The program I tried is posted below and you can test that the following command produces no output:

    python myCruzdbProgram.py | awk '/LCVSLC/'

How do I simply get the translation programatically for my region of interest?

    import cruzdb
    from cruzdb import Genome
    
    def translate_dna(sequence):
            gencode = { 'ATA':'I',
            'ATC':'I', 'ATT':'I', 'ATG':'M', 'ACA':'T', 'ACC':'T', 'ACG':'T',
            'ACT':'T', 'AAC':'N', 'AAT':'N', 'AAA':'K', 'AAG':'K', 'AGC':'S',
            'AGT':'S', 'AGA':'R', 'AGG':'R', 'CTA':'L', 'CTC':'L', 'CTG':'L',
            'CTT':'L', 'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCT':'P', 'CAC':'H',
            'CAT':'H', 'CAA':'Q', 'CAG':'Q', 'CGA':'R', 'CGC':'R', 'CGG':'R',
            'CGT':'R', 'GTA':'V', 'GTC':'V', 'GTG':'V', 'GTT':'V', 'GCA':'A',
            'GCC':'A', 'GCG':'A', 'GCT':'A', 'GAC':'D', 'GAT':'D', 'GAA':'E',
            'GAG':'E', 'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGT':'G', 'TCA':'S',
            'TCC':'S', 'TCG':'S', 'TCT':'S', 'TTC':'F', 'TTT':'F', 'TTA':'L',
            'TTG':'L', 'TAC':'Y', 'TAT':'Y', 'TAA':'\_', 'TAG':'\_', 'TGC':'C',
            'TGT':'C', 'TGA':'_', 'TGG':'W', }
            proteinseq = ''
            for n in range(0,len(sequence),3):
                    if gencode.has_key(sequence[n:n+3]) == True:
                            proteinseq += gencode[sequence[n:n+3]]
            return proteinseq
    
    g = cruzdb.Genome('hg38')
    
    # this is my region of interest
    chromosome = "chr19"
    start = 55014736
    stop = 55014753
    
    # retrieve annotation of the region with cruzdb
    gene = g.bin_query('refGene',chromosome,start,stop).all()
    
    # print the gene name ( it does find the correvt gene )
    print gene[0].name2
    
    # get the translated sequence
    sequence = "".join(gene[0].cds_sequence)
    print sequence
    
    # translation of my 18 base pair region should be LCVSLC but it is not
    print translate_dna(sequence.upper())
The first issue here is the coordinates you are interested in are on the negative strand, so you'll want to reverse complement the dna you obtain before trying to translate it.

Secondly, you can do all of this via UCSC tools like so:

    $ twoBitToFa -seq=chr19 -start=55014735 -end=55014753 http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit stdout  |  faRc stdin stdout | faTrans stdin stdout
    >RC_chr19:55014735-55014753
    CLSVCL

The twoBitToFa, faTrans and faRc commands are all available from our directory of utilities, in the directory appropriate for your operating system:<br>
http://hgdownload.soe.ucsc.edu/admin/exe/

If you are interested in querying the refGene table, it may be best to just [download it][1] and operate on it on locally your system so as not to violate our  Conditions of Use:<br>
http://genome.ucsc.edu/goldenPath/help/mysql.html

If you have further questions about UCSC data or tools feel free to send your question to one of the below mailing lists:

   - General questions: genome@soe.ucsc.edu
   - Questions involving private data: genome-www@soe.ucsc.edu
   - Questions involving mirror sites: genome-mirror@ose.ucsc.edu

ChrisL from the UCSC Genome Browser


  [1]: http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/refGene.txt.gz