Wednesday, August 12, 2009

Python simulation of Needleman-Wunsch 2

Continuing with the alignment problem, we have the BLOSUM50 matrix in a file titled BLOSUM.py (I can't find where I got it from, but there is a different version here).


#  Matrix made by matblas from blosum50.iij
# * column uses minimum score
# BLOSUM Clustered Scoring Matrix in 1/3 Bit Units
# Blocks Database = /data/blocks_5.0/blocks.dat
# Cluster Percentage: >= 50
# Entropy = 0.4808, Expected = -0.3573
A R N D C Q E G H I L K M F P S T W Y V B Z X *
A 5 -2 -1 -2 -1 -1 -1 0 -2 -1 -2 -1 -1 -3 -1 1 0 -3 -2 0 -2 -1 -1 -5
R -2 7 -1 -2 -4 1 0 -3 0 -4 -3 3 -2 -3 -3 -1 -1 -3 -1 -3 -1 0 -1 -5
N -1 -1 7 2 -2 0 0 0 1 -3 -4 0 -2 -4 -2 1 0 -4 -2 -3 4 0 -1 -5
D -2 -2 2 8 -4 0 2 -1 -1 -4 -4 -1 -4 -5 -1 0 -1 -5 -3 -4 5 1 -1 -5
C -1 -4 -2 -4 13 -3 -3 -3 -3 -2 -2 -3 -2 -2 -4 -1 -1 -5 -3 -1 -3 -3 -2 -5
Q -1 1 0 0 -3 7 2 -2 1 -3 -2 2 0 -4 -1 0 -1 -1 -1 -3 0 4 -1 -5
E -1 0 0 2 -3 2 6 -3 0 -4 -3 1 -2 -3 -1 -1 -1 -3 -2 -3 1 5 -1 -5
G 0 -3 0 -1 -3 -2 -3 8 -2 -4 -4 -2 -3 -4 -2 0 -2 -3 -3 -4 -1 -2 -2 -5
H -2 0 1 -1 -3 1 0 -2 10 -4 -3 0 -1 -1 -2 -1 -2 -3 2 -4 0 0 -1 -5
I -1 -4 -3 -4 -2 -3 -4 -4 -4 5 2 -3 2 0 -3 -3 -1 -3 -1 4 -4 -3 -1 -5
L -2 -3 -4 -4 -2 -2 -3 -4 -3 2 5 -3 3 1 -4 -3 -1 -2 -1 1 -4 -3 -1 -5
K -1 3 0 -1 -3 2 1 -2 0 -3 -3 6 -2 -4 -1 0 -1 -3 -2 -3 0 1 -1 -5
M -1 -2 -2 -4 -2 0 -2 -3 -1 2 3 -2 7 0 -3 -2 -1 -1 0 1 -3 -1 -1 -5
F -3 -3 -4 -5 -2 -4 -3 -4 -1 0 1 -4 0 8 -4 -3 -2 1 4 -1 -4 -4 -2 -5
P -1 -3 -2 -1 -4 -1 -1 -2 -2 -3 -4 -1 -3 -4 10 -1 -1 -4 -3 -3 -2 -1 -2 -5
S 1 -1 1 0 -1 0 -1 0 -1 -3 -3 0 -2 -3 -1 5 2 -4 -2 -2 0 0 -1 -5
T 0 -1 0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1 2 5 -3 -2 0 0 -1 0 -5
W -3 -3 -4 -5 -5 -1 -3 -3 -3 -3 -2 -3 -1 1 -4 -4 -3 15 2 -3 -5 -2 -3 -5
Y -2 -1 -2 -3 -3 -1 -2 -3 2 -1 -1 -2 0 4 -3 -2 -2 2 8 -1 -3 -2 -1 -5
V 0 -3 -3 -4 -1 -3 -3 -4 -4 4 1 -3 1 -1 -3 -2 0 -3 -1 5 -4 -3 -1 -5
B -2 -1 4 5 -3 0 1 -1 0 -4 -4 0 -3 -4 -2 0 0 -5 -3 -4 5 2 -1 -5
Z -1 0 0 1 -3 4 5 -2 0 -3 -3 1 -1 -4 -1 0 -1 -2 -2 -3 2 5 -1 -5
X -1 -1 -1 -1 -2 -1 -1 -2 -1 -1 -1 -1 -1 -2 -2 -1 0 -3 -1 -1 -1 -1 -1 -5
* -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 1



We load it with this code saved as BLOSUM.py:


def loadMatrix(fn='blosum50.txt'):
FH = open(fn,'r')
data = FH.read()
FH.close()
L = data.strip().split('\n')

# matrix has metadata lines beginning w/'#'
# also has extra rows,cols for 'BZX*'
L = [e for e in L if not e[0] in '#BZX*']

# the last 4 cols are also 'BZX*'
L = [e.split()[:-4] for e in L]
aaNames = L.pop(0)
# each row also starts with the AA name
L = [t[1:] for t in L]

M = dict()
for i in range(len(aaNames)):
for j in range(len(aaNames)):
k = aaNames[i] + aaNames[j]
M[k] = int(L[i][j])
return aaNames,M

def showM(M,aaNames):
#each element takes up 4 spaces
print ' ' + ' '.join(aaNames)
for i,aa1 in enumerate(aaNames):
print aa1,
L = list()
for j in range(i+1):
aa2 = aaNames[j]
k = aa1 + aa2
L.append(str(M[k]).rjust(4))
print ''.join(L)
for aa1 in aaNames:
for aa2 in aaNames:
assert M[aa1+aa2] == M[aa2+aa1],aa1+aa2

if __name__ == '__main__':
aaNames,M = loadMatrix(fn='blosum50.txt')
showM(M,aaNames)