Biopython: Enhance CodonAlignment to support addition with another alignment

[vanya-antonov]

Setup

I am reporting a problem with Biopython version, Python version, and operating
system as follows:

>>> import sys; print(sys.version)
3.5.2 (default, Nov 23 2017, 16:37:01)
[GCC 5.4.0 20160609]
>>> import platform; print(platform.python_implementation()); print(platform.platform())
CPython
Linux-4.4.0-135-generic-x86_64-with-Ubuntu-16.04-xenial
>>> import Bio; print(Bio.__version__)
1.71

Expected behaviour

Creation of a new MultipleSeqAlignment object by concatenating CodonAlignment with another CodonAlignment or MultipleSeqAlignment object via the + operator.

Actual behaviour

TypeError is produced in both cases

Steps to reproduce

from Bio.Alphabet import generic_dna, generic_protein
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Align import MultipleSeqAlignment
from Bio.codonalign import build
from Bio.codonalign.codonalphabet import default_codon_alphabet

# Generate a CodonAlignment object
seq1 = SeqRecord(Seq('ATGTCTCGT', alphabet=generic_dna), id='pro1')
seq2 = SeqRecord(Seq('ATGCGT', alphabet=generic_dna), id='pro2')
pro1 = SeqRecord(Seq('MSR', alphabet=generic_protein), id='pro1')
pro2 = SeqRecord(Seq('M-R', alphabet=generic_protein), id='pro2')
aln = MultipleSeqAlignment([pro1, pro2])
codon_aln = build(aln, [seq1, seq2])
print(codon_aln)

```
CodonAlphabet(Standard) CodonAlignment with 2 rows and 9 columns (3 codons)
ATGTCTCGT pro1
ATG---CGT pro2

```python
# Generate MultipleSeqAlignment object
tail1 = SeqRecord(Seq('AAA', alphabet=generic_dna), id='pro1')
tail2 = SeqRecord(Seq('AAA', alphabet=generic_dna), id='pro2')
tail_aln = MultipleSeqAlignment([tail1, tail2])
print(tail_aln)

DNAAlphabet() alignment with 2 rows and 3 columns
AAA pro1
AAA pro2

# Trying to concatenate them produces TypeError
# Should we allow concatenation of codon and non-codon alignments at all?
new_aln = codon_aln + tail_aln

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python3.5/dist-packages/Bio/Align/__init__.py", line 698, in __add__
    return MultipleSeqAlignment(merged, alpha, annotations, column_annotations)
  File "/usr/local/lib/python3.5/dist-packages/Bio/Align/__init__.py", line 177, in __init__
    self.extend(records)
  File "/usr/local/lib/python3.5/dist-packages/Bio/Align/__init__.py", line 525, in extend
    rec = next(records)
  File "/usr/local/lib/python3.5/dist-packages/Bio/Align/__init__.py", line 688, in <genexpr>
    merged = (left + right for left, right in zip(self, other))
  File "/usr/local/lib/python3.5/dist-packages/Bio/SeqRecord.py", line 859, in __add__
    answer = SeqRecord(self.seq + other.seq,
  File "/usr/local/lib/python3.5/dist-packages/Bio/Seq.py", line 306, in __add__
    return self.__class__(str(self) + str(other), a)
  File "/usr/local/lib/python3.5/dist-packages/Bio/codonalign/codonseq.py", line 76, in __init__
    raise TypeError("Input alphabet should be a CodonAlphabet object.")
TypeError: Input alphabet should be a CodonAlphabet object.

# Generate another CodonAlignment object
codon_tail_aln = MultipleSeqAlignment([tail1, tail2], alphabet=default_codon_alphabet)
print(codon_tail_aln)

CodonAlphabet(Standard) alignment with 2 rows and 3 columns
AAA pro1
AAA pro2

# Trying to concatenate them produces another TypeError
new_codon_aln = codon_aln + codon_tail_aln

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python3.5/dist-packages/Bio/Align/__init__.py", line 687, in __add__
    alpha = Alphabet._consensus_alphabet([self._alphabet, other._alphabet])
  File "/usr/local/lib/python3.5/dist-packages/Bio/Alphabet/__init__.py", line 396, in _consensus_alphabet
    alpha = Gapped(alpha, gap_char=gap)
  File "/usr/local/lib/python3.5/dist-packages/Bio/Alphabet/__init__.py", line 213, in __init__
    AlphabetEncoder.__init__(self, alphabet, gap_char)
  File "/usr/local/lib/python3.5/dist-packages/Bio/Alphabet/__init__.py", line 178, in __init__
    self.letters = alphabet.letters + new_letters
TypeError: can only concatenate list (not "str") to list

# A possible workaround is to convert CodonAlphabet to DNAAlphabet for
# all the sequences in the CodonAlignment object
for r in codon_aln._records:
    r.seq = Seq(str(r.seq), generic_dna)
new_aln = codon_aln + tail_aln
print(new_aln)

Gapped(Alphabet(), '-') alignment with 2 rows and 12 columns
ATGTCTCGTAAA pro1
ATG---CGTAAA pro2

[peterjc]

Supporting this looks like a very reasonable expectation, but it looks like it was never implemented/tested.

Here CodonAlignment inherits the addition method from MultipleSeqAlignment, but the codon specific alphabet rules in CodonAlphabet clash somewhere.

@zruan does this look like something you could work on?

[zruan]

Peter is correct. The codon alphabet clashes with the single letter alphabet rule.

In general, I don't think concatenating a CodonAlignment object with a MultipleSeqAlignment object is a good practice. If this is really necessary, converting the CodonAlignment object to a MultipleSeqAlignment object first is more reasonable. CodonAlignment class has a toMultipleSeqAlignment method for this purpose. The resulting concatenated MultipleSeqAlignment can then be converted back to CodonAlignment using a classmethod CodonAlignment.from_msa() with an appropriate codon alphabet.

Concatenating two CodonAlignment objects is reasonable. I can implement the function in the next a few days.

[peterjc]

Great. Do you think we can define:

• CodonAlignment + CodonAlignment -> CodonAlignment
• CodonAlignment + MultipleSeqAlignment -> MultipleSeqAlignment
• MultipleSeqAlignment + CodonAlignment -> MultipleSeqAlignment

[kaskales]

I see this issue is still open, is it ok if I give it a shot?

[zruan]

Hi @kaskales ,

Thank you for your interest. I implemented this function in my own fork of biopython, but I've never got a chance to write the test case. You can check out my code here. It will be good if you can write some unittest to the new add scheme.

Zheng

[kaskales]

@zruan Thank you for letting me know. I've pulled your changes and I will get started working on the unittests on my fork.