Python Programming for Biology: Bioinformatics and Beyond

Calculating sequence similarity

The  next  example  of  a  Python  function  will  consider  a  substitution  matrix  like  the  ones

discussed and use it to calculate an overall similarity score for two aligned sequences. The

inputs to the function are two strings of sequence letters and the similarity matrix.

def calcSeqSimilarity(seqA, seqB, simMatrix):

numPlaces = min(len(seqA), len(seqB))

totalScore = 0.0

for i in range(numPlaces):

residueA = seqA[i]

residueB = seqB[i]

totalScore += simMatrix[residueA][residueB]

---

return totalScore

# Test with pre-defined substitution matrices

# DNA example

print(calcSeqSimilarity('AGCATCGCTCT', 'AGCATCGTTTT', DNA_2))

# Protein example

print(calcSeqSimilarity('ALIGNMENT', 'AYIPNVENT', BLOSUM62))

The  calcSeqSimilarity()  function  is  very  similar  in  construction  to  the  previous

calcSeqIdentity() function, except that this time rather than seeing if two residue letters are

equal, we use them as keys to look up a similarity score in the substitution matrix. Note

that  this  function  has  one  big  deficiency:  it  cannot  deal  with  gaps  (‘-’).  To  address  this

problem  we  could  put  entries  for  gaps  into  the  similarity  matrix.  However,  a  simpler

solution  is  to  introduce  a  separate  gap  penalty;  gaps  are  generally  undesirable  but  are

tolerable if the subsequent alignment matches well. If we complicate things slightly more

we can have different gap penalties depending on whether we are inserting a new gap or

extending an existing one. Generally extending an existing gap (i.e. putting one dash after

another)  has  the  smaller  penalty.  This  is  equivalent  to  saying  that  we  score  alignments

more  highly  if  they  use  fewer,  longer  gapped  regions.  The  following  modified  function

uses gap penalties insert and extend, which both carry default values. Note that we have a

different  name  for  the  new  function  and  its  input  sequences  (alignA  and  alignB)  to

reinforce  the  fact  that  it  is  working  on  a  pair  or  aligned  sequences,  including  any  gaps,

rather than just plain sequences:

Pay  special  attention  to  the  logic  above  in  the  if/elif/else  statement.  If  a  gap  is  not

amongst the two residue codes then the score for the pair is obtained as before from the

similarity matrix dictionary. Otherwise, we do have a gap and thus carry on to check the

other two conditions. If the position in the sequence i is not at the very start (i > 0) and one

---

of the previous positions was a gap we subtract the extend penalty. And if all else fails we

have a gap that starts anew, so we subtract the insert penalty.