Global alignment with affinity

global_affinity_alignment.py

##########################
# Author : Pratik Anand  #
##########################
from tabulate import tabulate
import copy

class Obj :
	data = 0
	x = 0
	y = 0
	z = 0
        back="0"
	def __init__(self, data):
		self.data = data
   

MIN = -1000
USE_BLOSUM = False
match = 1
mismatch = -4
gap_open = -11
gap_extension = -1


def printTable( rowHeader, columnHeader, table ) :
	new_table =[]
	for innerRow in table :
		new_innerRow = [ ob.data for ob in innerRow ]
		new_table.append(new_innerRow)

	for indx,innerRow in enumerate(new_table):
		innerRow.insert( 0, rowHeader[indx])
	print tabulate( new_table, headers=columnHeader, tablefmt="grid" )



def sequence_align( strA, strB, charSet= None,blosum = None) :
	print strA + " vS " + strB
	strA = "0" + strA
	strB = "0" + strB

	megaL = []
	for i in range(len(strB)):
		innerL = [Obj(0) for x in range(len(strA))]   # list assignment evil
		megaL.append(innerL)

	megaX = []
	for i in range(len(strB)):
		innerX = [0 for x in range(len(strA))]   # list assignment evil
		megaX.append(innerX)

	megaY = []
	for i in range(len(strB)):
		innerY = [0 for x in range(len(strA))]   # list assignment evil
		megaY.append(innerY)


	for i in range( len(strB) ) :
		for j in range( len(strA) ) :
			
			if i==0 and j==0 :
				megaY[i][j] = 0
				megaX[i][j] = 0
				megaL[i][j].data = 0
	                        megaL[i][j].back="0"
				continue

			
			if (i==0 and j== 1) or ((i==1 and j== 0)) :
				megaY[i][j] = gap_open + gap_extension
				megaX[i][j] = gap_open + gap_extension
			else :
				#vertical
				megaX[i][j] = max(  ( megaX[i-1][j] if (i-1) >=0 else MIN ) + gap_extension,
									( megaL[i-1][j].data if (i-1) >=0 else MIN  ) + gap_open + gap_extension
								 )
				#horizantal
				megaY[i][j] = max(  ( megaY[i][j-1] if (j-1) >=0 else MIN ) + gap_extension,
									( megaL[i][j-1].data if (j-1) >=0 else MIN  ) + gap_open + gap_extension
								 )


			#Main TODO
			indx_i = -1
			indx_j = -1
			if USE_BLOSUM is True :
				listCharSet = list(charSet)
				for indx,string in enumerate(listCharSet) :
					if string == strB[i] :
						indx_i = indx
					if string == strA[j]:
						indx_j = indx

			local_matched = (match if strB[i] == strA[j] else mismatch) if USE_BLOSUM is False else (blosum[indx_i][indx_j] if blosum[indx_i][indx_j] != "?" else MIN )
			
			megaL[i][j].z = ( megaL[i-1][j-1].data if (i-1) >=0 and (j-1) >=0 else MIN ) + local_matched
			megaL[i][j].y = megaY[i][j]
			megaL[i][j].x = megaX[i][j]


			megaL[i][j].data = max( megaL[i][j].z, megaL[i][j].y, megaL[i][j].x )
	                if megaL[i][j].data == megaL[i][j].y :
	                        megaL[i][j].back="y"
	                elif megaL[i][j].data == megaL[i][j].x :
	                        megaL[i][j].back="x"
	                else : 
	                        megaL[i][j].back="z"







	#print-sandra
	for i in range( len(strB) ) :
		for j in range( len(strA) ) :
			print str(megaL[i][j].data) + " [Y:"+ str(megaL[i][j].y) + ",X:" + str(megaL[i][j].x)  + ",M:" + str(megaL[i][j].z) + "]  ",

		print "\n"
	

	printTable(strB, list(strA), megaL)

	ix=len(strB)-1
	iy=len(strA)-1
	seqA = []
	seqB = []
	while ix>0 and iy>0 :
	    if megaL[ix][iy].back == "z" :
	        seqA.append( strB[ix] )
	        seqB.append( strA[iy] )
	        ix = ix-1
	        iy = iy-1
	    elif megaL[ix][iy].back == "y" :
	        seqA.append("_" )
	        seqB.append( strA[iy] )
	        iy = iy-1
	    else :
	        seqB.append("_")
	        seqA.append( strB[ix] )
	        ix = ix-1
	        
	seqA.reverse()
	seqB.reverse()
	print "".join(seqA)
	print "".join(seqB)


if __name__ == '__main__' :

	s1 = "GATG"
	s2 = "GATATC"

	sequence_align(s1, s2)