It works! Thank you very much!
• 0 views
•
link
written
by
Charles Yin
How to display mismatched sequences from alignment when using Biopython
When using Biopython alignment to align two protein sequences, the output does not show the mismatched amino acids. Would you please have some suggestions to display the mismatched amino acids by Biopython?
The following is the code for alignment by Biopython. I tried different options, the problem could not be resolved.
----------------------------------------------------------
from Bio import pairwise2
from Bio.pairwise2 import format_alignment
proteinSeqQuery='HRTLGLLLHTQVSIQQLLKLPAECFHPKPKVNSVLIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKHAK'
proteinSeqseqFound='HRTLGLLLHTQVSIKQLLKLPAECFHPKPKVNSALIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKYAK'
alignments = pairwise2.align.localxx(proteinSeqQuery,proteinSeqseqFound)
print('Score',alignments[0][2])
print(format_alignment(*alignments[0]))
----------------------------------------------------------------
The program output and the output from EMBI-EBI EMBOSS Water (Protein Alignment) are shown as follows:
Biopython:
HRTLGLLLHTQVSI-QQLLKLPAECFHPKPKVNSVLIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKHAK
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
HRTLGLLLHTQVSIKQ-LLKLPAECFHPKPKVNSALIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKYAK
EMBL-EBI:
EMBOSS_001 1 HRTLGLLLHTQVSIQQLLKLPAECFHPKPKVNSVLIKLTRHTTDVPDKYW 50
||||||||||||||:||||||||||||||||||.||||||||||||||||
EMBOSS_001 1 HRTLGLLLHTQVSIKQLLKLPAECFHPKPKVNSALIKLTRHTTDVPDKYW 50
EMBOSS_001 51 KLYTYFVSKWVNREYRQLFTKNQFHQAMKHAK 82
|||||||||||||||||||||||||||||:||
EMBOSS_001 51 KLYTYFVSKWVNREYRQLFTKNQFHQAMKYAK 82
and in this picture: output
Thank you in advance!
• 11,197 views
•
link
updated
by
st.ph.n
•
written
by
Charles Yin
•
written
by
Charles Yin
2 answers
Here's a bit of a bodge that might work. Note though that it doesn't take in to consideration amino acid similarity. It would be easy enough to code in with a few more elif statements though provided you can put together a table of what matches with what.
In this instance it just leaves a blank where there's no match, by simply iterating both strings from pairwise2 and checking equivalency. You could create a slightly fancier print output too if desired but I think this works.
In essence all you need to do is create a custom string printed in the middle of the 2 alignments. It's generated by appending to a list, and then joining the list to make another string. Standard use of monospaced fonts in terminals etc means that the alignments of all the character columns is taken care of for you.
from Bio import pairwise2
from Bio.pairwise2 import format_alignment
proteinSeqQuery='HRTLGLLLHTQVSIQQLLKLPAECFHPKPKVNSVLIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKHAK'
proteinSeqseqFound='HRTLGLLLHTQVSIKQLLKLPAECFHPKPKVNSALIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKYAK'
alignments = pairwise2.align.localxx(proteinSeqQuery,proteinSeqseqFound)
match = []
for a, b in zip(alignments[0][0],alignments[0][1]):
if a == b:
match.append('|')
else:
match.append(' ')
print(alignments[0][0])
print("".join(match))
print(alignments[0][1])
Yields:
HRTLGLLLHTQVSIQ-QLLKLPAECFHPKPKVNSV-LIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKH-AK
|||||||||||||| |||||||||||||||||| ||||||||||||||||||||||||||||||||||||||||||||| ||
HRTLGLLLHTQVSI-KQLLKLPAECFHPKPKVNS-ALIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMK-YAK
Also I'd maybe play around with the gap penalty and mismatch penalty a little bit since that alignment seems to be scoring gaps quite low.
For instance, it looks to me like the 2 gaps here should instead just be a mismatch, but I could be wrong...
QVSIQ-QLLKL
|||| |||||
QVSI-KQLLKL
Would make more sense to me as:
QVSIQQLLKL
|||| |||||
QVSIKQLLKL
Thank you! yes, that is a problem. What gap penalty do you use in localxx() for the correct result (second one)? It looks localms() function should be used. The following is the improved version.
from Bio import pairwise2
seqQuery='HRTLGLLLHTQVSIQQLLKLPAECFHPKPKVNSVLIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKHAK'
seqFound='HRTLGLLLHTQVSIKQLLKLPAECFHPKPKVNSALIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKYAK'
# Four parameters (2,-1,-1,-1) in alignment:Identical characters are given 2 points, 1 point is deducted for each non-identical character.
# 1 point is deducted when opening a gap, and 1 point is deducted when extending it.
alignments = pairwise2.align.localms(seqQuery,seqFound, 2, -1, -1, -1)
match = []
for a, b in zip(alignments[0][0],alignments[0][1]):
if a == b:
match.append('|')
else:
match.append('*')
m="".join(match)
s=[]
s.append(alignments[0][0]+'\n')
s.append(m+'\n')
s.append(alignments[0][1])
alignedSeqs="".join(s)
print('\n')
print(alignedSeqs)
The output is as follows:
HRTLGLLLHTQVSIQQLLKLPAECFHPKPKVNSVLIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKHAK
||||||||||||||*||||||||||||||||||*|||||||||||||||||||||||||||||||||||||||||||||*||
HRTLGLLLHTQVSIKQLLKLPAECFHPKPKVNSALIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKYAK
• 0 views
•
link
written
by
Charles Yin
Thanks! There are a lot of confusions in the alignment functions in Biopython. I currently use Biopython v1.67, and will upgrade to v1.68.
• 0 views
•
link
written
by
Charles Yin
Have a look in Biopython's API. pairwise2 is complex and allows many options. To get the same results as in EMBOSS' WATER you should use the same parameters, e.g. using align.localds(seq1, seq2, matrix, gap-open-penalty, gap-extend-penalty). In the API is an example how to implement the matrix.
Is there a way to make this faster? It seems really slow right now.
• 0 views
•
link
written
by
biogirl04 •
#!/usr/bin/env python
query='HRTLGLLLHTQVSIQQLLKLPAECFHPKPKVNSVLIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKHAK'
target='HRTLGLLLHTQVSIKQLLKLPAECFHPKPKVNSALIKLTRHTTDVPDKYWKLYTYFVSKWVNREYRQLFTKNQFHQAMKYAK'
aln = ''
for n in range(len(query)):
if query[n] == target[n]:
aln += '.'
else:
aln += query[n]
print '{}\t{}'.format('Query', query)
print '{}\t{}'.format('Aln', aln)
print '{}\t{}'.format('Target', target)
Thanks for this try, but it may not handle gaps in a sequence alignment.
• 0 views
•
link
written
by
Charles Yin
Log in to answer this question.
I added (code) markup to your post for increased readability. You can do this by selecting the text and clicking the 101010 button. When you compose or edit a post that button is in your toolbar, see image below:
Thanks! I will use this approach for formatting code.
Please say which version of Biopython you are using, since the pairwise code was rewritten in Biopython 1.68