Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # returns the cost of a substitution of characters a -> b
- def substitute(a, b):
- aIndex = letterList.index(a)
- bIndex = letterList.index(b)
- blosumVal = blosumMatrix[aIndex][bIndex]
- return blosumVal
- # return the alignment and direction matrix for Needleman-Wunsch algorithm
- def nwAlignment(stringA, stringB, gapPenalty):
- # get matrix dimensions
- m = len(stringA) + 1
- n = len(stringB) + 1
- # initalise the alignment and direction matrix
- alignMatrix = np.zeros((n,m), dtype=int)
- directionMatrix = np.zeros((n,m), dtype=str)
- # insert the border values (all gaps vertical and horizontal)
- for i in range(n):
- alignMatrix[i][0] = gapPenalty * i
- directionMatrix[i][0] = 'V'
- for j in range(m):
- alignMatrix[0][j] = gapPenalty * j
- directionMatrix[0][j] = 'H'
- # overwrite the corner value in direction matrix
- directionMatrix[0][0] = '-'
- # calculate the remaining matrix values as the max of the 3 backtracking options
- for i in range(1, n):
- for j in range(1, m):
- a = stringA[j-1]
- b = stringB[i-1]
- vGap = alignMatrix[i-1][j] + gapPenalty
- hGap = alignMatrix[i][j-1] + gapPenalty
- sub = alignMatrix[i-1][j-1] + substitute(a, b)
- bestMove = max(sub, hGap, vGap)
- alignMatrix[i][j] = bestMove
- # insert the character in the direction matrix to identify which option was the best move
- if bestMove == sub:
- directionMatrix[i][j] = 'D'
- elif bestMove == hGap:
- directionMatrix[i][j] = 'H'
- else:
- directionMatrix[i][j] = 'V'
- return (alignMatrix, directionMatrix)
- # returns the 2 global matching strings
- def getGlobalMatch(stringA, stringB, directionMatrix):
- # inialises the start position as the bottom right corner of the matrix
- i = len(stringB)
- j = len(stringA)
- # recursive function that builds up the match strings
- def nextMove(i, j, outputA, outputB):
- direction = directionMatrix[i][j]
- if direction == '-':
- outputA = outputA[::-1]
- outputB = outputB[::-1]
- return (outputA, outputB)
- elif direction == 'D':
- outputA = outputA + stringA[j-1]
- outputB = outputB + stringB[i-1]
- result = nextMove(i-1, j-1, outputA, outputB)
- return result
- elif direction == 'H':
- outputA = outputA + stringA[j-1]
- outputB = outputB + '-'
- result = nextMove(i, j-1, outputA, outputB)
- return result
- else:
- outputA = outputA + '-'
- outputB = outputB + stringB[i-1]
- result = nextMove(i-1, j, outputA, outputB)
- return result
- result = nextMove(i, j, '', '')
- return result
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement