I have pairs of coding DNA sequences which I wish to perform pairwise codon alignments via Python, I have "half completed" the process.

So far..

  • I retrive pairs of orthologous DNA sequences from genbank using Biopython package.
  • I translate the orthologous pairs into peptide sequences and then align them using EMBOSS Needle program.

I wish to..

  • Transfer the gaps from the peptide sequences into the original DNA sequences.


I would appreciate suggestions for programs/code (called from Python) that can transfer gaps from aligned peptide sequence pairs onto codons of the corresponding nucleotide sequence pairs. Or programs/code that can carry out the pairwise codon alignment from scratch.

enter image description here

  • $\begingroup$ Why don't you simply align the nucleotide sequences in the first place? Just use any alignment software. $\endgroup$ – terdon Dec 30 '13 at 17:59
  • $\begingroup$ Codon alignments serve a purpose, only this way you can obtain values of non-synonymous (dN) and synonymous (dS) rates of mutation.. which can then be used to obtain the dN/dS ratio used to infer selection $\endgroup$ – hello_there_andy Dec 30 '13 at 18:04
  • $\begingroup$ You can do that easily enough from the DNA sequence, no need to translate. Just align and read the codons from the alignment. You seem to be complicating your life needlessly. That said, this is probably too computer geeky for this site and you would be better off asking on biostars. $\endgroup$ – terdon Dec 30 '13 at 18:07
  • 5
    $\begingroup$ This question appears to be off-topic because it is about programming and not biology. Bioinformatics questions are on topic here if they are about the underlying concepts, you are asking a very technical issue whose answer will be code. $\endgroup$ – terdon Dec 30 '13 at 18:10
  • 1
    $\begingroup$ I am quite familiar with it, thanks :). I did my PhD on comparative genomics and reconstructing evolutionary histories, I've done more alignments than I care to remember. You need to know whether two codons are synonymous or not in order to get the dN/dS but not necessarily to translate them, it's just an easy way to check. Anyway, using my dict suggestion this is pretty trivial. Feel free to @ping me in Biology Chat if you need more help with it. $\endgroup$ – terdon Dec 30 '13 at 18:30

The basic process would be (in pseudocode, I don't know python well enough, I'm a Perl geek):

for ($i=0;$i<length($seq1);$i++){
for ($i=0;$i<length($seq2);$i++){

At this point you will have two arrays or hashes or tuples or dicts or whatever holding the amino acid that corresponds to each codon position of the input seqs. You go and do your translation and dS/dN thing and you can then use these lists to map back to the original nucleotides:

for ($i=0;$i<length($ProteinSseq1);$i++){
  print codons1[$i]

You will need to adapt this to deal with codons correctly (i++ will only increment by one, you need to extract the triplets) and then to read the gaps but these are all programming issues that have nothing to do with the biological side of things.

| improve this answer | |

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