Question: 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())

Answer: 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:

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:

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