Needleman-Wunsch Algorithm In Ruby

• Tweet
• 

        As described at http://en.wikipedia.org/wiki/Needleman-Wuncsh

# uses MDArray
require 'mdarray'

def needle(sequence, reference)
	gap = -5
	
	# similarity matrix
	#
	#     A  G  C  T
	#  A 10 -1 -3 -4
	#  G -1  7 -5 -3
	#  C -3 -5  9  0
	#  T -4 -3  0  8
	
	s = { 'AA' => 10,
	        'AG' => -1,
	        'AC' => -3,
	        'AT' => -4,
	        'GA' => -1,
	        'GG' =>  7,
	        'GC' => -5,
	        'GT' => -3,
	        'CA' => -3,
	        'CG' => -5,
	        'CC' =>  9,
	        'CT' =>  0,
	        'TA' => -4,
	        'TG' => -3,
	        'TC' =>  0,
	        'TT' =>  8 }
	
	rows = reference.length + 1
	cols = sequence.length + 1
	a = MDArray.new(rows, cols)
	for i in 0...(rows) do a[i,0] = 0 end
	for j in 0...(cols) do a[0,j] = 0 end
	for i in 1...(rows)
		for j in 1...(cols)
			choice1 = a[i-1, j-1] + s[(reference[i-1].chr + sequence[j-1].chr).upcase]
			choice2 = a[i-1, j] + gap
			choice3 = a[i, j-1] + gap
			a[i,j] = [choice1, choice2, choice3].max
		end
	end
			
	ref = ''
	seq = ''
	i = reference.length 
	j = sequence.length
	while (i > 0 and j > 0)
		score = a[i,j]
		score_diag = a[i-1, j-1]
		score_up = a[i, j-1]
		score_left = a[i-1, j]
		if (score == score_diag + s[reference[i-1].chr + sequence[j-1].chr])
			ref = reference[i-1].chr + ref
			seq = sequence[j-1].chr + seq
			i -= 1
			j -= 1
		elsif (score == score_left + gap)
			ref = reference[i-1].chr + ref
			seq = '-' + seq
			i -= 1
		elsif (score == score_up + gap)
			ref = '-' + ref
			seq = sequence[j-1].chr + seq
			j -= 1
		end
	end
	
	while (i > 0)
		ref = reference[i-1].chr + ref
		seq = '-' + seq
		i -= 1		
	end
	
	while (j > 0)
		ref = '-' + ref
		seq = sequence[j-1].chr + seq
		j -= 1
	end
	
	[seq, ref] 
end

Tags:
• alignment
• Bioinformatics
• Ruby