Counting letters in phylip alignment columns with Biopython

Question

I have been using python 3.6 and biopython 1.72 to work with protein data files. I am using a protein sequence file (phylip format), for example:

    14    678
Zebrafish  LSSCGVVSGD LISVILPASS LEETQTSSAA AHQTHTDQQA GGSHVSSSSS 
Fugu       LASCGIVSGD LIRVALIRAA LMAADAPDRA DDGGGHSEQA VAMATNQVS- 
Chicken    LASCGIVSGD LIRVALIRAA LMAADAPDRA DDGGGHSEQA VAMATNQVS- 
Zebra      LASYGIVPGD LICLLLEETD AKSSLPPPSS SSPPSLQNGR ETSTLTPKNQ

I want to manipulate the file in a way so that;

1. I may get the number of columns which have all the same letters (amino acids)
2. The number of columns which have at least 2 different amino acids (with no dashes).
3. The number of columns which have any single dash.

I have iterated over the columns of the file with:

alignment = AlignIO.read(open("example.phy"), "phylip") 
def column(alignment):
    n = len(alignment[0].seq)

    for i in range(n):
        c = []

        for record in alignment:
            c.append(record.seq[i])

        yield c

This returns the columns of the file, but I don't understand how I can manipulate them further.

Can I get some helpful suggestions?

Answer 1

I suggest you use zip() to transform your data into columns and then apply the following counts similar to heathobrien suggests:

from Bio import AlignIO

alignment = AlignIO.read(open('example.phy'), 'phylip')

var1 = var2 = var3 = 0
for col in zip(*alignment):
    num_unique = len(set(col))
    if num_unique == 1 and col[0] != '-':
        var1 += 1
    elif num_unique > 1 and '-' not in col:
        var2 += 1
    elif '-' in col:  # assumes 1 or more dashes
        var3 += 1

print('Number of columns with the same amino acid: {}\n'
      'Number of columns with at least 2 amino acids (no gaps): {}\n'
      'Number of columns with at least one gap: {}'
       .format(var1, var2, var3))

For your example data, this outputs:

Number of columns with the same amino acid: 9
Number of columns with at least 2 amino acids (no gaps): 40
Number of columns with at least one gap: 1

Answer 2

using something like for col in column(alignment):

1. I may get the number of columns which have all the same letters (amino acids)

   len(set(col)) == 1
   

2. The number of columns which have at least 2 different amino acids (with no dashes):

   len(set(col)) > 1 and '-' not in col
   

3. The number of columns which have any single dash:

   '-' in col
   

or if you really want columns with a single dash:

    col.count('-') == 1