3
$\begingroup$

I need to convert contigs into their respective protein sequences given a reference genome (i.e. I need to take a substring, whose position is already known on the string, and I need to locate the nearest start and stop codons).

This is tricky because sometimes the first position of triplet codon within the contig will not be a multiple of three and could be located in the intergenic region (i.e. between genes). The goal is to separate both coding and non-coding DNA.

This is my code so far:

from Bio.Seq import Seq
from Bio.Alphabet import generic_dna, generic_protein

start_codons = ['ATG']
stop_codons = ['TAG', 'TAA', 'TGA']

string = 'GG*TAG*CCAATT*ATG*AACGAA*TAG*GAC' #remove '*', just for visual
contigs = ['CCAA', 'TGAAC', 'GAA', 'GGAC']
positions = [5, 12, 17, 22] #position indices for each contig on string

extended_contigs = []
extended_position_contigs = []
intergenic_contigs = []
intergenic_position_contigs = []
for i in contigs:
    extended_contigs.append(#some code)
    extended_positions_contigs.append(#some code)
    intergenic_contigs.append(#some code)
    intergenic_positions_contigs.append(#some code)

I should get extended_contigs = ['ATGAAC', 'ATGAACGAA'] and extended_positions_contigs = [12, 17]. These are the contigs that are located within genes. In order to code them into peptides, I need to reach back in the string until I find the start codon and extend the initial contig (e.g. TGAAC -> ATGAAC and GAA -> ATGAACGAA)

I should also get intergenic_contigs = ['CCAA', 'GGAC'] and intergenic_positions_contigs = [5, 22]. When the missing code is run, the computer searches to the left of the string and finds a stop codon (e.g. TAG) before the start codon. Thus, the contig is located in between two genes and nothing needs to be added. These intergenic contigs are just stored in a new list.

My code continues:

prot_contigs = []
for i in extended_contigs:
    my_dna = Seq(i, generic_dna)
    my_prot = my_dna.translate()
    prot_contigs.append(str(my_prot))

Here, no new code needs to be added. After the above is run, prot_contigs = ['MN', 'MNE'].

The last step of code (which I need help with) will transform prot_contigs into new_prot_contigs = ['MN', 'E'].

How? If for any contig (e.g. 'TGAAC'), the beginning or end is part of another codon (not perfect multiple of 3), the extra codon on either end will be maintained (e.g. 'MN' stays 'MN'). Otherwise, if the contig (e.g. 'GAA') perfectly matches a codon, anything that was added onto it will be removed (e.g. 'MNE' becomes 'E').

I would try to solve the 2 portions of code myself, but I am unsure on how to take a position on a string (i.e. beginning points of contig) and look along the string to find the nearest start/stop codons, so I could determine the DNA's function and accurately sequence the protein coding contigs into peptides.

Any help would be appreciated!

Improve this question
$\endgroup$
3
  • $\begingroup$ Are you trying to predict open reading frames and you don’t want to use someone else’s gene prediction tool? $\endgroup$
    – Chris_Rands
    Oct 19, 2018 at 14:36
  • $\begingroup$ @Chris_Rands Unfortunately, I am not familiar with open reading frames or gene prediction tools. Would such a software be able to solve the above? If so, is there open source code I could use? $\endgroup$
    – Thomas Ribeiro
    Oct 19, 2018 at 17:48
  • 1
    $\begingroup$ @ThomasRibeiro you can explore the code of TransDecoder which is pretty widely used. $\endgroup$
    – NatWH
    Oct 21, 2018 at 12:51

1 Answer 1

Reset to default
1
$\begingroup$

If your string is 'GGTAGCCAATTATGAACGAATAGGAC' then I would assume that you could separate this into coding and non-coding regions such as:

GG*TAG*CCAATT*    (non-coding)
ATG*AACGAA*TAG    (coding region)
*GAC              (non-coding)

It starts with a start codon ATG and ends with an end codon TAG.

Your first event position 12 does not make sense to me. What does * mean? It seems of that your sequence starts with GG which is not a triplet? Perhaps did you mean GGG?

I'm the author of node based DNA visualisation AminoSee. It has a bug in that it only renders half right: because DNA can be read backwards, and my program only does it forwards. I actually came to BioInformatics today to try and get my head around reading backwards and how that might work.

AminoSee only transcodes 3 base pairs at a time, and discards triplets that aren't, well, 3 base pairs longs. I'm self taught and so I guess unaware if it's possible to have any sequence that is not a multiple of 3. Perhaps it's possible. I assumed not. My program is more of an art piece that a science tool, but I plan to use it one day to perform a "visual diff" between two species (Humans and Chimp). See AminoSee DNA Visualisation over at aminosee.funk.nz

Improve this answer
$\endgroup$

Your Answer

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

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