If 2 protein sequence X of length m and Y of length n and if there are several highest scoring global alignments, then I want to get one alignment that has largest number of columns in which a letter from X is aligned to an identical letter from Y. I need an algorithm to find alignment with maximum number of identical matches. Can anyone help me?
1 Answer
There are a lot of possible ways to count the number of perfect matches in an alignment, and the details will depend on which alignment tool you are using.
But, if you use the built-in pairwise aligner in BioPython (Bio.Align.PairwiseAligner()
, https://biopython.org/docs/1.75/api/Bio.Align.html), then the easiest way is to use the formatted string representation of the alignment, which uses the pipe/vertical-bar character ('|') to indicate perfect matches.
Example code:
from Bio import Align
aligner = Align.PairwiseAligner()
alignments = aligner.align("GAACT", "GAT")
max_matches = 0
for i, a in enumerate(alignments):
align_str = f'{a}'
align_str = align_str.split('\n')
match_str = align_str[1]
match_count = match_str.count('|')
mismatch_count = match_str.count('.')
print(f'Number of matches: {match_count}')
print(f'Number of mismatches: {mismatch_count}')
if match_count>max_matches:
max_matches = match_count
max_ind = i
best_alignment = alignments[max_ind]
In this example, both alignments have the same number of perfect patches (3), but that may not always be the case, depending on the input sequences and the scoring parameters used. The example just uses the default scoring parameters, but in general, you probably want to choose parameters appropriate for your type of alignment.