#!/usr/bin/ruby -w ################################################################################# # File : gap.rb # # Author : Jean-Eudes Duchesne (duchesnj@iro.umontreal.ca) # # Created : 2 janvier 2004 # # Last updated : 3 février 2004 # # Description : Implementation of the global alignment algorithm for comparing # # two DNA sequences using a ponderation matrix. # # # # Usage : global.rb S T # # S et T are the strings to be compared. The script is pretty # # basic, so try not to use extreme cases. For exemple, output # # formatting is pretty simple, using strings that are too long # # might give unexpected results ! Also, only works for sequences # # in lower case. # # # ################################################################################# ### Libs ### ### Globals ### Ws = -2 # penalty for starting a gap We = -0.1 # penalty for extending a gap ### Classes & Methods ### class Array2D def initialize @data = [[]] end def [](i,j) # Want i to be the rows and j the columns if @data[j] == nil || @data[j][i] == nil return nil else return @data[j][i] end end def []=(i,j,x) @data[j] = [] if @data[j] == nil @data[j][i] = x end end class PondMatrix def initialize # Not very O.O., but I need the values ! @data = Array2D.new @data[0,0] = @data[1,1] = @data[2,2] = @data[3,3] = 3 @data[0,1] = @data[1,0] = @data[0,2] = @data[2,0] = @data[1,3] = @data[3,1] = @data[2,3] = @data[3,2] = -2 @data[0,3] = @data[3,0] = @data[1,2] = @data[2,1] = -1 end def [](a,b) #printf("A-%s-A\n",a) #printf("B-%s-B\n",b) case a when "a" i = 0 when "c" i = 1 when "g" i = 2 when "t" i = 3 end case b when "a" j = 0 when "c" j = 1 when "g" j = 2 when "t" j = 3 end return @data[i,j] end end def MAX(x,y) unless(y == nil || x == nil) if(y > x) return y else return x end end if(y == nil) return x else return y end end ### Main ### S = ARGV[0] T = ARGV[1] m = S.length + 1 n = T.length + 1 # Initialisation of ponderation matrix p = PondMatrix.new # Initialisation of arrays D = Array2D.new E = Array2D.new F = Array2D.new for i in 0 ... m D[i,0] = F[i,0] = Ws + (i*We) end for j in 0 ... n D[0,j] = E[0,j] = Ws + (j*We) end D[0,0] = E[0,0] = F[0,0] = 0 # calculate every other cell for i in 1 ... m for j in 1 ... n E[i,j] = MAX(D[i-1,j]+Ws, E[i-1,j]) + We F[i,j] = MAX(D[i,j-1]+Ws, F[i,j-1]) + We D[i,j] = MAX(D[i-1,j-1],MAX(E[i-1,j-1], F[i-1,j-1])) + p[S[i-1..i-1],T[j-1..j-1]] end end #Print dynamic tables and minimal nb of operations printf("\n--- D[i,j] ---") printf("\n - ") for j in 0 ... n printf(" %s ", T[j..j]) end printf("\n- ") for i in 0 ... m for j in 0 ... n printf("%5.1f ", D[i,j]) end printf("\n%s ", S[i..i]) end printf("\n--- E[i,j] ---") printf("\n - ") for j in 0 ... n printf(" %s ", T[j..j]) end printf("\n- ") for i in 0 ... m for j in 0 ... n printf("%5.1f ", E[i,j]) end printf("\n%s ", S[i..i]) end printf("\n--- F[i,j] ---") printf("\n - ") for j in 0 ... n printf(" %s ", T[j..j]) end printf("\n- ") for i in 0 ... m for j in 0 ... n printf("%5.1f ", F[i,j]) end printf("\n%s ", S[i..i]) end printf("\nScore : %6.1f\n\n", D[m-1,n-1])