High-throughput classification of amino acid differences in protein sequence alignment

Question

[aerobiceukaryote]

Hi,

I am looking for a way to take a pairwise alignment and output the amino acid modification. For example, the sequences could be:

>seq1
MYPROT

>seq2
MYYROT

The desired output would be: P3Y

I have hundreds of protein sequences so I'd like something high-throughput because doing this by eye is tedious and error-prone.

Any help is greatly appreciated!

Answer 1

[Mensur Dlakic]

This may not be exactly what you want, but it should give you a solid starting point. Below are two sequences (maybe save them into test.fas), of which the second is mutated randomly in 7 positions compared to first:

>rseq_000000
CLDHEASRFAMPMQYDRAVEQCCNIHVNIDTLPASRAQLGQYTDPTGGVF
LPHPDRVAEGDCRQCNLDDVKSTGVMEFSCNGQSTSPSSSAVTSPVHACL
VMGSSFHYNEAALRTHNLDHEEFCATAITDSSQDQKIAKTELVQESYRED
ALIGMKEPRAQPVGRPRDPCASPVKRRKLVPNIVSTHYKCACQVKKNCPD
EPTRNGLFRVEGLNEGKPGSQIVGGGRSKGVQSSIIDSPTSSEGGIPRVV
PPGLSPSPLETISREECHAQETLIDRRSLQQAEYNAPHQYPMIKEPPGVQ
YDMKLQEELQTDRACFDQPSQDGYPMPDSLSSPGEAV
>rseq_000000_mut
CLDHNASRFAMPMQYDRAVEQCCNIHVNIDTLPASRAQLGQYTDPTGGVF
LPHSDRVAEGDCEQCNVDDVKSTGVMEFSCNGQSTSPSSSAVTSPVHACL
VMGSSFHYNEAALRTHNLDHEEFCATAITDSSQDQKIAKTELVQESYRED
ALIGMKEPRAQPVGRPREPCASPVKRRKLVPNIVSTHYKCACQVKKNCPD
EPTRNALFRVEGLNEGKPGSQIVGGGRSKGVQSSIIDSPTSSEGGIPRVV
PPGLSPSPIETISREECHAQETLIDRRSLQQAEYNAPHQYPMIKEPPGVQ
YDMKLQEELQTDRACFDQPSQDGYPMPDSLSSPGEAV

This short python code may be saved into file fasta_diff.py:

import sys
import numpy as np
from Bio import SeqIO

FastaFile = open(sys.argv[1], 'r')
store = np.array([])

for seqs in SeqIO.parse(FastaFile, 'fasta'):
    store = np.append(store, str(seqs.seq))

for z in range(len(store[0])):
    if store[0][z] != store[1][z]:
        print('%s' % (store[0][z]+str(z+1)+store[1][z]))

FastaFile.close()

Running python fasta_diff.py test.fas will produce the following output:

E5N
P54S
R63E
L67V
D168E
G206A
L259I

[aerobiceukaryote]

Thank you Mensur!