score_pairwise.py

"""
Code written to answer stackoverflow question on substitution matrices:
http://stackoverflow.com/questions/5686211/
"""

from Bio.SubsMat import MatrixInfo

def score_match(pair, matrix):
    """ 
    Return score for a given pair of residues in a give matrix.
    
    Matrices only have one triangle, so if res1 -> res2 isn't in matrix 
    res2 -> res1 will be... 
    """
    if pair not in matrix:
        return matrix[(tuple(reversed(pair)))]
    else:
        return matrix[pair]

def score_pairwise(seq1, seq2, matrix, gap_s, gap_e):
    """Score a pairwise alignment, with gap_s and gap_e as the gap penalties """
    score = 0
    gap = False
    for i in range(len(seq1)):
        pair = (seq1[i], seq2[i])
        if not gap:
            if '-' in pair:
                gap = True
                score += gap_s
            else:
                score += score_match(pair, matrix)
        else:
            if '-' not in pair:
                gap = False
                score += score_match(pair, matrix)
            else:
                score += gap_e
    return score

def test():
  """Test these functions on the provided alignment and the blosum 62 matrix"""
  seq1 = 'PAVKDLGAEG-ASDKGT--SHVVY----------TI-QLASTFE'
  seq2 = 'PAVEDLGATG-ANDKGT--LYNIYARNTEGHPRSTV-QLGSTFE'
  score =  score_pairwise(seq1, seq2, MatrixInfo.blosum62, -5, -1)
  assert score == 82
  print "score: %s" % score

if __name__ == "__main__":
  test()