Question: Need help with my Needleman-Wunsch Python program (traceback) (error)
0
gravatar for daniel.merlinx
6 months ago by
daniel.merlinx0 wrote:

Good Afternoon Biostar Community,

I have a problem with my python program. I tried to implement the Needleman-Wunsch algorithm but I get the wrong alignment due to a mistake in the traceback section I assume. I got the correct matrix but not the correct sequence alignment. With this NW("SCHULFACH","SCHLUESSEL","2","-1","-1") input values I should get this solution "SCH_ULFACH || SCHLUESSEL" but I get this "----SCHULFACH and SCHL---UESSEL". It would help me a lot if you have a look at my program and try to figure out what I've done wrong.

Thanks in advance!

def NW(Seq1,Seq2,match,mismatch,gap):

#Initialization of the matrix
matrix=[[x for x in range(len(Seq1)+1)] for i in range(len(Seq2)+1)]

#set matrix to zero
for x in range(len(Seq1)+1):
    for i in range(len(Seq2)+1):
        matrix[i][x]=0



#Calculation of the matrix
traceback=[[x for x in range(len(Seq1)+1)] for i in range(len(Seq2)+1)]
a=0
b=0
if len(Seq1)==len(Seq2):
    a=len(Seq1)+1
    b=len(Seq2)+1
else:
    if len(Seq1)>len(Seq2):
        a=len(Seq2)+1
        b=len(Seq1)+1

    else:
        b=len(Seq1)+1
        a=len(Seq2)+1





for i in range(a):
    for x in range(b):

     if x==0 and i==0:
        matrix[0][0]=0

     else:

        if x==0:
            f1=float("-inf")
            f2=matrix[i-1][x]
            f3=float("-inf")

        elif i==0:
            f1=float("-inf")
            f2=float("-inf")
            f3=matrix[i][x-1]

        else:
            f1=matrix[i-1][x-1]
            f2=matrix[i][x-1]
            f3=matrix[i-1][x]



        if(x>0 and i>0):
            if(Seq2[i-1]==Seq1[x-1]):
                m=int(match)
            else:
                m=int(mismatch)
        else:
            m=int(mismatch)



        w1=f1+m
        w2=f2+int(gap)
        w3=f3+int(gap)
        erg=(w1,w2,w3)
        w4=max(erg)


        if w4==w1:
            traceback[i][x]=0

        elif w4==w2:
            traceback[i][x]=1 


        else:#w4==w3
            traceback[i][x]=-1


        try:
            matrix[i][x]=w4

        except IndexError:
            print("ERROR")

#traceback

alnseqA=""#horizontal
alnseqB=""#vertical
j=len(Seq1)#horizontal
i=len(Seq2)#vertical
print(j,i)

while(i>0 or j>0):
    if traceback[i][j]==0:
        alnseqA += Seq1[j-1]
        alnseqB += Seq2[i-1]
        j=j-1
        i=i-1
    elif traceback[i][j]==1:
        alnseqA += "-"
        alnseqB += Seq2[i-1]
        i=i-1
    else:
        alnseqB += "-"
        alnseqA += Seq1[j-1]
        j=j-1
print(matrix)
print(alnseqA[::-1])
print(alnseqB[::-1])


return (alnseqA,alnseqB)

NW("GCTCAA","CTAAG","2","-1","-2")

NW("SCHULFACH","SCHLUESSEL","2","-1","-1")

ADD COMMENTlink written 6 months ago by daniel.merlinx0

It seems you're using the wrong scores:

if x==0:
            f1=float("-inf")
            **f2=matrix[i-1][x]** 
            f3=float("-inf")

        elif i==0:
            f1=float("-inf")
            f2=float("-inf")
            **f3=matrix[i][x-1]**

        else:
            f1=matrix[i-1][x-1]
            **f2=matrix[i][x-1]** # should be f2=matrix[i-1][x]
            **f3=matrix[i-1][x]** # should be f3=matrix[i][x-1]

edit: fixed mis-annotation

ADD REPLYlink modified 6 months ago • written 6 months ago by cschu1811.9k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1271 users visited in the last hour