How to iterate protein sequences using amino acids?

Question

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
    Canis           REWSSARPERSKGRRKPVDAAAVSAVQTSQTSSDVAVSSSCRSMEMQDLTSPHSRLSGSS
    Mus             -------------------------------------------MEMQDLTSPHSRLSGSS
    Rattus          -------------------------------------------MEMQDLTSPHSRLSGSS
                                                               *****************

Example phylip file:

     3   120
Apis       MPWIKCLGGG GGKARRGKPL SVSQPIHLLV KSEFEPQCHV FRTKQLRKPR 
Branchiost ---------- ---------- ---------- ---------- ---------- 
Nematostel ---------- ---------- ---------- ---------- ---------- 

           PCHLCHQAVI KQASCCRVCK YICHKTCEDK GARFRDTSQQ IHKQWQADPA 
           ---------- ---------- ---------- ---------- ---------- 
           ---------- ---------- ---------- ---------- ---------- 

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?

Answer 1

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).