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 contig will not be a multiple of three (i.e. the first three nucleotides in the contig may not perfectly match a codon). Also, sometimes the contigs 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()

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!

  • $\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 '18 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 '18 at 17:48
  • $\begingroup$ @ThomasRibeiro you can explore the code of TransDecoder which is pretty widely used. $\endgroup$ – NatWH Oct 21 '18 at 12:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.