rendered paste body#! /usr/bin/pythonfrom numpy import *import sysimport argparse #Requires Python 3.2 or the attached argparse.py fileclass NW_Align: def __init__(self,str1,str2): self.match = "|" self.mismatch = " " self.gap = "-" self.matchScore = 2 self.mismatchScore = -1 self.gapScore = -2 x = len(str1) + 1 y = len(str2) + 1 self._str1 = str1.lower() self._str2 = str2.lower() self._matrix = zeros(x*y).reshape(x,y) for i in range(x): if i > 0: self._matrix[i,0] = self._matrix[i-1,0] - 2 for i in range(y): if i > 0: self._matrix[0,i] = self._matrix[0,i-1] - 2 for i in range(len(str1)): for j in range(len(str2)): self._alignPair(i+1,j+1) def calcAlignments(self, max): score = self._matrix[len(self._str1),len(self._str2)] paths = self._calcPaths(max) alignmentCount = len(paths) output = str(alignmentCount) + " Optimal Alignment" if alignmentCount > 1: output += "s" output += " found:\n\n" for i in range(alignmentCount): alignmentList = self._toString(paths[i]) output += "Alignment No. " + str(i + 1) + "\n" for j in range(3): for site in alignmentList: output += site[j] output += "\n" output += "\n" return output def _printMatrix(self): print self._matrix def _alignPair(self,x,y): a = self._str1[x-1] b = self._str2[y-1] pairwise = 2 if (a == b) else -1 pairwise += self._matrix[x-1,y-1] topGap = self._matrix[x-1,y] - 2 leftGap = self._matrix[x,y-1] - 2 results = array([pairwise, topGap, leftGap]) max = results.max() self._matrix[x,y] = max def _toString(self, path): alignment = [] for i in range(len(path)-1): site = () dX = path[i+1][0] - path[i][0] dY = path[i+1][1] - path[i][1] X = path[i][0] Y = path[i][1] if dX == 0 and dY == 1: site = (self.gap,self.mismatch,self._str2[Y]) elif dX == 1 and dY == 0: site = (self._str1[X],self.mismatch,self.gap) elif dX == 1 and dY == 1: if self._str1[X] == self._str2[Y]: site = (self._str1[X],self.match,self._str2[Y]) else: site = (self._str1[X],self.mismatch,self._str2[Y]) alignment.append(site) return alignment def _calcPaths(self, max): path = [(len(self._str1),len(self._str2))] queue = [path] complete = [] if max > 0: while (len(queue) > 0): current = queue.pop(0) size = len(queue) + len(complete) newSteps = self._findSteps(current, size, max) if len(newSteps) > 0: for step in newSteps: newPath = list(current) newPath.insert(0, step) queue.append(newPath) else: complete.append(current) return complete def _findSteps(self, path, size, max): options = [] lastStep = path[0] lastX = lastStep[0] lastY = lastStep[1] if size + len(options) < max and self._isMisOrMatch(lastX, lastY): options.append((lastX-1,lastY-1)) if size + len(options) < max and self._isTopGap(lastX, lastY): options.append((lastX-1,lastY)) if size + len(options) < max and self._isLeftGap(lastX, lastY): options.append((lastX,lastY-1)) return options def _isMisOrMatch(self,x,y): result = False if x > 0 and y > 0: stepCost = self._matrix[x,y] - self._matrix[x-1,y-1] if self._str1[x-1] == self._str2[y-1]: result = (stepCost == self.matchScore) else: result = (stepCost == self.mismatchScore) return result def _isTopGap(self,x,y): result = False if x > 0: stepCost = self._matrix[x,y] - self._matrix[x-1,y] result = (stepCost == self.gapScore) return result def _isLeftGap(self,x,y): result = False if y > 0: stepCost = self._matrix[x,y] - self._matrix[x,y-1] result = (stepCost == self.gapScore) return result if __name__ == "__main__": parser = argparse.ArgumentParser(description='This is a Needleman-Wunsch alignment program') parser.add_argument("-m", "--max", action='store', dest='max', type=int, default=1, help='Maximum number of alignments to output, default value: 1') parser.add_argument("--first", action='store', dest='first', type=str, default=None, help='First Sequence (Required)', required=True) parser.add_argument("--second", action='store', dest='second', type=str, default=None, help='Second Sequence (Required)', required=True) options = parser.parse_args() array = NW_Align(options.first, options.second) results = array.calcAlignments(options.max) print results