Python simulation of Needleman-Wunsch 5

To finish up with alignments, the code I posted previously can be called with a command-line argument that is the name of the target file containing two FASTA-formatted sequences.

Here is one of my favorite comparisons---I use it often in class. We compare the protein product of the yeast CDC28 gene with the human CDK2 gene. I always mention that the human gene complements the yeast cell cycle defect, and comment that "the cell cycle was obviously invented a long time ago."

Here is the output:

MSGELANYKRLEKVGEGTYGVVYKALDLRPGQGQRVVALKKIRLESEDEG M E N++++EK+GEGTYGVVYKA + + G+ VVALKKIRL++E EG M--E--NFQKVEKIGEGTYGVVYKARN-K-LTGE-VVALKKIRLDTETEG VPSTAIREISLLKELKDDNIVRLYDIVHSDAHKLYLVFEFLDLDLKRYME VPSTAIREISLLKEL NIV+L D++H++ +KLYLVFEFL DLK++M+ VPSTAIREISLLKELNHPNIVKLLDVIHTE-NKLYLVFEFLHQDLKKFMD GIP-KDQPLGADIVKKFMMQLCKGIAYCHSHRILHRDLKPQNLLINKDGN PL ++K ++ QL +G+A+CHSHR+LHRDLKPQNLLIN +G ASALTGIPL--PLIKSYLFQLLQGLAFCHSHRVLHRDLKPQNLLINTEGA LKLGDFGLARAFGVPLRAYTHEIVTLWYRAPEVLLGGKQYSTGVDTWSIG +KL DFGLARAFGVP+R YTHE+VTLWYRAPE+LLG K YST VD WS+G IKLADFGLARAFGVPVRTYTHEVVTLWYRAPEILLGCKYYSTAVDIWSLG CIFAEMCNRKPIFSGDSEIDQIFKIFRVLGTPNEAIWPDIVYLPDFKPSF CIFAEM R+ +F GDSEIDQ+F+IFR LGTP+E +WP + +PD+KPSF CIFAEMVTRRALFPGDSEIDQLFRIFRTLGTPDEVVWPGVTSMPDYKPSF PQWRRKDLSQVVPSLDPRGIDLLDKLLAYDPINRISARRAAIHPYFQE-S P+W R+D+S+VVP LD G LL ++L YDP RISA+ A HP+FQ+ + PKWARQDFSKVVPPLDEDGRSLLSQMLHYDPNKRISAKAALAHPFFQDVT

Python simulation of Needleman-Wunsch 4

Finally, here is the code to organize everything. If run without a command-line argument, it does three short comparisons. Here is the output:

TGCTCGTA T TC TA T--TCATA HEAGAWGHE-E AW HE E --P-AW-HEAE HKKLYLVFEFLDLD-RYMEGIPKE + ++YLVFE+LDL+ ++M+ P++ N-RIYLVFEYLDLETKFMDSCPED

And here is the code:

import string,sys from BLOSUM import loadMatrix,showM from Utils import pprint,report from Utils import printAlignment,load from Inits import newDict,setUp,bases from Inits import seqType,initScore from Algorithm import doScoring,trackback def run(s1,s2,kind=None,v=False): if not kind: kind = seqType(s1,s2) sc = initScore() # set up the box L = setUp(s1,s2,sc) if kind == 'protein': # m = blosum aaNames,m = loadMatrix( fn='blosum50.txt') #showM(m,aaNames) else: m = dict() for n1 in bases: for n2 in bases: k = n1 + n2 if n1 == n2: m[k] = 5 else: m[k] = -2 doScoring(L,s1,s2,m,sc) t = trackback(L,s1,s2,m) if not v: print '\n'.join(t) + '\n' else: printAlignment(t) def main(): s1 = 'TGCTCGTA' s2 = 'TTCATA' run(s1,s2) s1 = 'HEAGAWGHEE' s2 = 'PAWHEAE' run(s1,s2) s1 = 'HKKLYLVFEFLDLDRYMEGIPKE' s2 = 'NRIYLVFEYLDLETKFMDSCPED' run(s1,s2) if __name__ == "__main__": fn = None try: fn = sys.argv[1] except IndexError: pass if fn: s1,s2 = load(fn) run(s1,s2,v=True) else: main()

Python simulation of Needleman-Wunsch 3

This is the third of four posts showing code for a simulation to do Needleman-Wunsch alignments. In this installment, we get the algorithm itself. It is obvious by inspection that the doScoring function implements the NW algorithm. Check out the function handlePos in the trackback function.

Save this code as Algorithm.py in the same location as the other files.

from Inits import newDict def doScoring(L,s1,s2,matrix,sc): R = len(s1) + 1 # row length C = len(s2) + 1 # col length # for each column for c in range(2, R+1): # for each row in that column for r in range(2, C+1): i = (r-1)*R + c-1 up = L[i-R] if up['path'] == 'D': upscore = up['score'] + sc.gap else: upscore = up['score'] + sc.ext left = L[i-1] if left['path'] == 'D': leftscore = left['score'] + sc.gap else: leftscore = left['score'] + sc.ext diag = L[i-R-1]['score'] m = s1[c-2] n = s2[r-2] diag += matrix[m+n] # for debugging # report(r,c,i,upscore,leftscore,diag) if (diag >= leftscore) and (diag >= upscore): L[i] = newDict(diag, 'D') elif (leftscore > upscore): L[i] = newDict(leftscore, 'L') else: L[i] = newDict(upscore, 'U') def trackback(L,s1,s2,blosum): R = len(s1) + 1 # items per row or numCols def handlePos(i,s1L,s2L): j,k = i%R-1,i/R-1 D = L[i] #print D['score'],i,j,k,s1[j],s2[k] if D['path'] == 'U': s1L.append('-') s2L.append(s2[k]) return i-R if D['path'] == 'L': s1L.append(s1[j]) s2L.append('-') return i-1 if D['path'] == 'D': s1L.append(s1[j]) s2L.append(s2[k]) return i-(R+1) s1L = list() s2L = list() i = len(L) - 1 while i > 0: i = handlePos(i,s1L,s2L) s1L.reverse() s2L.reverse() mL = list() for i,c1 in enumerate(s1L): c2 = s2L[i] if '-' in c1 or '-' in c2: mL.append(' ') elif c1 == c2: mL.append(c1) elif blosum[c1+c2] > 0: mL.append('+') else: mL.append(' ') retL = [''.join(s1L),''.join(mL),''.join(s2L)] return retL

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)

Python simulation of Needleman-Wunsch 1

I want to show you my Python code for doing Needleman-Wunsch global alignments. I hasten to remind you that I don't really expect anyone to study this stuff, but I want to encourage you to try the problem yourself. It's fun, and will increase your understanding of the method. The purpose of posting the code is to give you a working example to fall back on if you encounter difficulty.

Let's do some (boring) utility functions first. I have the following code in a file called Inits.py. The setUp function makes a matrix of the type we used in the demo:

bases = 'ACGT' def initScore(g=-6,e=-6): class Score: pass score = Score() score.gap = g score.ext = e return score def seqType(s1,s2): for c in s1+s2: if not c in bases: return 'protein' return 'DNA' def newDict(score, path): return {'score':score, 'path':path} def setUp(s1,s2,sc): R = len(s1) + 1 # items per row or numCols C = len(s2) + 1 # items per col or numRows # a list of D with keys = score,path # just use a flat list, not array L = [None]*R*C L[0] = { 'score':0,'path':None } L[1] = newDict(sc.gap, 'L') for c in range(2,R): score = L[c-1]['score'] + sc.ext L[c] = newDict(score, 'L') L[R] = newDict(sc.gap, 'U') for r in range(2, C): prev = (r-1)*R next = r*R score = L[prev]['score'] + sc.ext L[next] = newDict(score, 'U') #print len(L) return L

And I have functions to load the sequences (filename can be passed in on the command line), and for pretty printing of alignments, saved in the following module as Utils.py:

# utility: pretty print alignment def pprint(L,s1,s2): N = len(s1)+1 # row length pad = 5 # print title row using s1 print ' '*pad*2, for c in s1: print c.rjust(pad), print print i = 0 L = L[:] while L: # print col 1 using chars from s2 if i == 0: print ' '.rjust(pad), else: print s2[i-1].rjust(pad), i += 1 sub = L[:N] L = L[N:] for D in sub: if D: score = str(D['score']) path = D['path'] if path: print (score+path).rjust(pad), else: print score.rjust(pad), else: print '.'.rjust(pad), print print #=============================================== # for debugging def report(r,c,i,up,left,diag): print 'r', r, 'c' ,c, print 'i',i if up: print 'up', up if left: print 'left', left if diag: print 'diag', diag print s1[c-2] print s2[r-2] def printAlignment(t): print L1,L2,L3 = t SZ = 50 for i in range(0,len(L1),SZ): print L1[i:i+SZ] print L2[i:i+SZ] print L3[i:i+SZ] print def load(fn): FH = open(fn,'r') data = FH.read() FH.close() if not data.strip()[0] == '>': return 'data is not FASTA' sys.exit() try: s1,s2 = data.strip().split('\n\n') except ValueError: return 'need two sequences to compare' sys.exit() s1 = s1.split('\n',1)[1] s2 = s2.split('\n',1)[1] s1 = s1.replace('\n','') s2 = s2.replace('\n','') return s1,s2

Alignment: Gotoh, and Smith-Waterman

In sequence alignment, we want to penalize gaps heavily so that we don't get this:



However, insertions and deletions (indels) don't happen one nucleotide at a time. Instead, it is thought that they tend to occur in the surface exposed loops of proteins, rather than in the core. Core regions are tightly packed and unlikely to tolerate much disruption. I need to look for papers about the above assertions, but let's continue. The idea of affine gap penalties is that there is a large hit for allowing a gap in the alignment, but a rather small penalty for extending one. We might have:

• gap opening penalty = -12
• gap extension penalty = -2

In a real problem, we would follow this approach. However, I want to keep it simple, so for this post, we are going to do as before, and have one penalty for each character (base or amino acid) opposite a gap (-8).

We want to align protein sequences this time: HEAGAWGHEE and PAWHEAE. We'll use one of the BLOSUM scoring matrices to do this (Henikoff & Henikoff, PMID: 1438297. Here, I have copied out the relevant pairwise scores from the BLOSUM50 matrix:



We set up the table for the alignment scores exactly as before:



We fill it in using the same rules:



Having cached the arrows which show how we got to each square, we can do the traceback starting in the lower right-hand corner and work our way back:



The alignment is then:



As before, notice that going up corresponds to a gap in sequence 1, while going left is a gap in sequence 2.

This particular application of NW is said to be due to Gotoh (PMID: 7166760). His paper in JMB is buried in the bowels of the central library about 2 miles from my desk. So I haven't actually read it. If you have a copy, I'd be grateful.

Smith-Waterman is a modification of this approach. It is an algorithm for local sequence alignment. It is the same as before, but with a simple new idea: if the accumulated score goes negative, set it equal to zero. And start the traceback from the maximum score:



This optimization eliminates the noise of poorly matched segments. SSearch (here) is a commonly used implementation. It is still computationally expensive (O(n²). We will need something faster.

Alignment: Needleman-Wunsch

Sergey called and said that if I don't do any bioinformatics examples they are going to make me change the name of the blog. :)

Here is a first installment.

The Needleman-Wunsch algorithm (PMID: 5420325---if this reference doesn't make any sense to you, just enter that number into the search box on this page), carries out a global alignment on two sequences. The wikipedia article on this topic is not very clear, in my opinion. The NW method produces a global rather than a local alignment. It is guaranteed to be the best possible (highest scoring) alignment. It is not guaranteed to be biologically relevant.

I found a nice explanation in Dave Swofford's lectures. Maybe you can find them, but his employer FSU has taken them down, apparently. That's too bad, because he helped me understand lots of things, most important, the methods for Hidden Markov Models, like the Viterbi algorithm, etc., which I hope to get to in a future post.

[In a perfect world, we would all put our lecture notes on line. But there are problems: our employers seem to think they "own" these presentations. And then there are people like me, who shamelessly use graphics snatched off the web to illustrate our slides. Due to stupid copyright issues, I can't put them up for you to see.]

There is a nice page about NW here. The example I will use is from Sean Eddy (PMID 15229554). He has several other very nice tutorials as well.

Let's say we want to align these two sequences: TGCTCGTA and TTCATA. We need scoring rules to evaluate each position in the alignment. Suppose we take:

• match = +5
• mismatch = -2
• insertion = -6

Notice this is an empirical procedure. We choose parameters which make the ultimate solutions seem correct. Set up a table like this:



At each turn, we move either down, right, or both right and down (on the diagonal). We start with an extra column at the left and an extra row at the top. Fill out these extra columns as shown. Moving right corresponds to an alignment with a gap in the sequence TTCATA, while moving down corresponds to a accepting a gap in the other sequence. Each additional base opposite a gap accumulates a penalty of -6 in our scoring system.

We calculate a running score for each position. In the top row, there is only one way to get to a square (from the left). Hence we add an additional -6 to the score for each move across the top. Similar logic applies for the left column.

Having finished the preliminaries, the first interesting position is the square where the initial nucleotide of each sequence lines up. We make the following calculation. We can do one of three things:

• move down, add -6 to the score from the square above
• move right, add -6 to the score from the square to the left
• move diagonally, and take the score from the upper left
in this case, evaluate the match or mismatch to find the score



We calculate these values, and take the maximum value. The mathematicians have a fancy way of writing that:



(These graphics are from Durbin et al, another very fine book).

We write the result of our calculation into the square, and also write down an arrow that shows which path we took.



In this case, we got +5 for the match, while the potential moves down or right had scores of -12. We continue to fill in the table. Here, I show the result, but I have not shown all of the arrows.



The last step is to begin at the right hand bottom corner, and follow the arrows in reverse back to the beginning. This is called the traceback. The resulting alignment is:



Notice the correspondence of the gaps in the top sequence with horizontal arrows in the table.

NW is a dynamic programming algorithm. We:

• break a big problem into smaller problems
• these have solutions that can be built up in stages
• keep track of (memoize) solutions to the small problems

• gives an alignment which is guaranteed to be the "best"
• will happily align unrelated sequences
• is computationally expensive: O(n²)