Bio.codonalign package

Submodules

Module contents

Code for dealing with Codon Alignments.

Bio.codonalign.build(pro_align, nucl_seqs, corr_dict=None, gap_char='-', unknown='X', codon_table=None, complete_protein=False, anchor_len=10, max_score=10)

Build a codon alignment from protein alignment and corresponding nucleotides.

Arguments:
  • pro_align - a protein MultipleSeqAlignment object

  • nucl_seqs - an object returned by SeqIO.parse or SeqIO.index or a collection of SeqRecord.

  • corr_dict - a dict that maps protein id to nucleotide id

  • complete_protein - whether the sequence begins with a start codon

Return a CodonAlignment object.

The example below answers this Biostars question: https://www.biostars.org/p/89741/

>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> from Bio.Align import MultipleSeqAlignment
>>> from Bio.codonalign import build
>>> seq1 = SeqRecord(Seq('ATGTCTCGT'), id='pro1')
>>> seq2 = SeqRecord(Seq('ATGCGT'), id='pro2')
>>> pro1 = SeqRecord(Seq('MSR'), id='pro1')
>>> pro2 = SeqRecord(Seq('M-R'), id='pro2')
>>> aln = MultipleSeqAlignment([pro1, pro2])
>>> codon_aln = build(aln, [seq1, seq2])
>>> print(codon_aln)
CodonAlignment with 2 rows and 9 columns (3 codons)
ATGTCTCGT pro1
ATG---CGT pro2

Using the newer codon aligner in Bio.Align, this analysis can be performed as follows (note that here we are also calculating the amino acid alignment between the two protein sequences):

>>> from Bio.Align import global_align, codon_align
>>> seq1 = SeqRecord(Seq('ATGTCTCGT'), id='pro1')
>>> seq2 = SeqRecord(Seq('ATGCGT'), id='pro2')
>>> pro1 = SeqRecord(Seq('MSR'), id='pro1')
>>> pro2 = SeqRecord(Seq('MR'), id='pro2')
>>> protein_alignment = global_align(pro1, pro2)[0]
>>> print(protein_alignment)
pro1              0 MSR 3
                  0 |-| 3
pro2              0 M-R 2

>>> alignment1 = codon_align(pro1, seq1)[0]
>>> alignment2 = codon_align(pro2, seq2)[0]
>>> print(alignment1)
pro1              0 M  S  R   3
pro1              0 ATGTCTCGT 9

>>> print(alignment2)
pro2              0 M  R   2
pro2              0 ATGCGT 6

>>> codon_alignment = protein_alignment.mapall([alignment1, alignment2])
>>> print(codon_alignment)
pro1              0 ATGTCTCGT 9
                  0 |||---||| 9
pro2              0 ATG---CGT 6

For comparison, this is what you get if you directly align the nucleotide sequences to each other:

>>> naive_alignments = global_align(seq1, seq2)
>>> for naive_alignment in naive_alignments:
...     print(naive_alignment)
pro1              0 ATGTCTCGT 9
                  0 |||-|--|| 9
pro2              0 ATG-C--GT 6

pro1              0 ATGTCTCGT 9
                  0 |||---||| 9
pro2              0 ATG---CGT 6

i.e. the correct alignment but also an incorrect alignment with incomplete codons.