The Needleman-Wunsch algorithm is a way to align sequences in a way that optimizes "similarity". Usually, a grid is generated and then you follow a path down the grid (based off the largest value) to compute the optimal alignment between two sequences. I have created a Python program, that given two strings, will create the resulting matrix for the Needleman-Wunsch algorithm. Currently, the program when ran will generate two random sequences of DNA and then print out the resulting Needleman-Wunsch matrix.
needlemanwunsch.py
import random
from tabulate import tabulate
class NeedlemanWunsch:
"""
Class used for generating Needleman-Wunsch matrices.
"""
def _compute_block(self, result, i, j):
"""
Given a block (corresponding to a 2 x 2 matrix), calculate the value o-
f the bottom right corner.
(Based on the equation:
M_{i,j} = max(M_{i-1,j-1} + S_{i,j},
M_{i,j-1} + W,
M_{i-1,j} + W)
Found here: https://vlab.amrita.edu/?sub=3&brch=274&sim=1431&cnt=1)
Args:
result : The current matrix that is being computed.
i : The right most part of the block being computed.
j : The bottom most part of the block being computed.
Returns:
The value for the right bottom corner of a particular block.
"""
return max(result[i-1][j-1] +
self._calc_weight(self._second_seq[i-1],
self._first_seq[j-1]),
result[i-1][j] + self.gap,
result[i][j-1] + self.gap)
def _calc_weight(self, first_char, second_char):
"""
Helper function, given two characters determines (based on the sc-
oring scheme) what the score for the particular characters can be.
Args:
first_char : A character to compare.
second_char : A character to compare.
Returns:
Either self.match or self.mismatch.
"""
if first_char == second_char:
return self.match
else:
return self.mismatch
def generate(self, first_seq, second_seq):
"""
Generates a matrix corresponding to the scores to the Needleman-Wu-
nsch algorithm.
Args:
first_seq : One of the sequences to be compared for similarity.
second_seq : One of the sequences to be compared for
similarity.
Returns:
A 2D list corresponding to the resulting matrix of the Needlem-
an-Wunsch algorithm.
"""
# Internally requies that the first sequence is longer.
if len(second_seq) > len(first_seq):
first_seq, second_seq = second_seq, first_seq
self._first_seq = first_seq
self._second_seq = second_seq
# Adjust sequence with "intial space"
# Initialize the resulting matrix with the initial row.
result = [list(range(0, -len(first_seq) - 1, -1))]
# Create initial columns.
for i in range(-1, -len(second_seq) - 1, -1):
row = [i]
row.extend([0]*len(first_seq))
result.append(row)
# Sweep through and compute each new cell row-wise.
for i in range(1, len(result)):
for j in range(1, len(result[0])):
result[i][j] = self._compute_block(result, i, j)
# Format for prettier printing.
for index, letter in enumerate(second_seq):
result[index + 1].insert(0, letter)
result[0].insert(0, ' ')
result.insert(0, list(" " + first_seq))
return result
def __init__(self, match=1, mismatch=-1, gap=-1):
"""
Initialize the Needleman-Wunsch class so that it provides weights for
match (default 1), mismatch (default -1), and gap (default -1).
"""
self.match = match
self.mismatch = mismatch
self.gap = gap
self._first_seq = ""
self._second_seq = ""
def deletion(seq, pos):
"""
Deletes a random base pair from a sequence at a specified position.
Args:
seq : Sequence to perform deletion on.
pos : Location of deletion.
Returns:
seq with character removed at pos.
"""
return seq[:pos] + seq[pos:]
def base_change(seq, pos):
"""
Changes a random base pair to another base pair at a specified position.
Args:
seq : Sequence to perform base change on.
pos : Locaion of base change.
Returns:
seq with character changed at pos.
"""
new_base = random.choice("ACTG".replace(seq[pos], ""))
return seq[:pos] + new_base + seq[pos:]
def mutate(seq, rounds=3):
"""
Mutates a piece of DNA by randomly applying a deletion or base change
Args:
seq : The sequence to be mutated.
rounds : Defaults to 3, the number of mutations to be made.
Returns:
A mutated sequence.
"""
mutations = (deletion, base_change)
for _ in range(rounds):
pos = random.randrange(len(seq))
seq = random.choice(mutations)(seq, pos)
return seq
def main():
"""
Creates a random couple of strings and creates the corresponding Needleman
-Wunsch matrix associated with them.
"""
needleman_wunsch = NeedlemanWunsch()
first_seq = ''.join(random.choices("ACTG", k=5))
second_seq = mutate(first_seq)
data = needleman_wunsch.generate(first_seq, second_seq)
print(tabulate(data, headers="firstrow"))
if __name__ == '__main__':
main()
I ended up using a NeedlemanWunsch
class, because using only function resulted in a lot DRY for the parameters match
, mismatch
, and gap
.
I am not particularly fond of the code. I haven't used numpy
or any related libraries because I couldn't see any way that it would significantly shorten the code, however, I would be willing to use numpy
if there is a significantly shorter way of expressing the matrix generation. However, the code seems frail, and very prone to off by one errors.