I was working on certain program using python and I have been using a protein aligned sequence file in two formats, phylip (.phy) and clustal (.aln). Example clustal file:

CLUSTAL 2.1 multiple sequence alignment

    Homo            -------------------------------------------MEMQDLTSPHSRLSGSS
    Macaca          -------------------------------------------MEMQDLTSPHSRLSGSS
    Mus             -------------------------------------------MEMQDLTSPHSRLSGSS
    Rattus          -------------------------------------------MEMQDLTSPHSRLSGSS

Example phylip file:

     3   120
Branchiost ---------- ---------- ---------- ---------- ---------- 
Nematostel ---------- ---------- ---------- ---------- ---------- 

           ---------- ---------- ---------- ---------- ---------- 
           ---------- ---------- ---------- ---------- ---------- 

I want to navigate through the file in a way that ;

  1. If the columns of the file have same amino acid, it saves and prints the amino acid that exists in the column.
  2. I need to know how can I search the columns by using the amino acids, if particular group of amino acids exist in that column or not. For example, if 'STA'or 'HY' or 'FVLIM' or 'NDEQHK' (with no '-') exists in any of the columns, then put a '@' in the file at the end of that column. Following is the code I have been trying to manipulate the file with the following program:

Code below

alignment = open ("example.aln")
a = b = c = 0 
for col in (alignment): 
    num = len(set(col)) 
        if num == 1 and col[0] != '-':
            print (num)
        a += 1
        elif num > 1 and '-' not in col:
        print (num)
        b += 1 elif '-' in col:  # assumes 1 or more dashes
        c += 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 one/more gaps: {}'
        .format(a, b, c))

These variables return the occurrences of amino acids in the columns of the file but I don't understand how to search using the amino acids. What is the best way to iterate over the file and write into a new file with results using python? Should I use python switch statement?

  • $\begingroup$ Hi Sidra, as it was the block of code was unreadable since it seemed like many of the lines had been merged together. Could you edit your post to make sure that the code displayed above exactly matches your program? It would also help if you make your code a working minimal example where you provide a small testcase "example.aln". meta.stackexchange.com/questions/22186/… $\endgroup$
    – conchoecia
    Sep 8 '18 at 15:08
  • $\begingroup$ Cross posted at Biostars: biostars.org/p/336237/#337089 $\endgroup$
    – Joe Healey
    Sep 10 '18 at 18:46

In response to question 1, save the following as a python3 script (e.g., aln_parser.py):

#!/usr/bin/env python3    

from Bio import AlignIO
import sys

aln = sys.argv[1]

algnmnt = AlignIO.read(open(aln, 'r'), 'clustal')

for col in range(0, algnmnt.get_alignment_length()):
    proteins = set(algnmnt[:,col])
    if len(proteins) == 1:
        print(f"Postion {col}: protein {''.join(proteins)}")

Run the script to query your alignment (clustal.aln).

python3 aln_parser.py clustal.aln

To save the output to file:

python3 aln_parser.py clustal.aln > results.txt

Requires python3 and Biopython (pip3 install biopython).

  • 1
    $\begingroup$ Looks good. You could also add a hash-bang #!/bin/env python3 to make it executable. $\endgroup$
    – Kamil S Jaron
    Oct 2 '18 at 15:17
  • $\begingroup$ question 2 has still been the problem, if anyone can help.! $\endgroup$ Oct 4 '18 at 7:52
  • $\begingroup$ @SidraYounas I would try to help, but I don't fully understand the question. $\endgroup$ Oct 26 '18 at 4:10
  • $\begingroup$ @user3479780 You can see the code sample above, that I have used. In this, I am iterating over the columns of the files. All I was trying to use the way if I can search for the presence of some particular amino acids for each column. $\endgroup$ Oct 26 '18 at 5:42
  • $\begingroup$ @user3479780 I have been trying to delete these columns which have a single amino acid, but I am not being successful. Can you help? $\endgroup$ Dec 3 '18 at 8:21

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.