5
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ BTW alignment and sequence-alignment tags should probably be synonyms $\endgroup$ – Chris_Rands Aug 30 '18 at 9:38
3
$\begingroup$

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
$\endgroup$
4
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Earlier I have used a for loop to iterate over the columns for c in column(alignment): and get the columns with dashes if "-" in c: list2.append(c) and without dashes else: list3.append(c) in two separate lists. But when it comes to count the total number of columns with dashes, it results out the number of dashes for each of the column in file and keeps on printing unless all the columns are iterated (because of for loop). $\endgroup$ – Sidra Younas Aug 30 '18 at 5:26
  • $\begingroup$ It would be helpful if you can elaborate me what you have suggested because I cannot follow up the idea. I still have my problem. 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. – $\endgroup$ – Sidra Younas Aug 30 '18 at 6:41

Your Answer

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

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