#!/usr/bin/ruby -w ################################################################################# # File : local.rb # # Author : Jean-Eudes Duchesne (duchesnj@iro.umontreal.ca) # # Created : 2 février 2004 # # Last updated : 2 février 2004 # # Description : Implementation of the local alignment algorithm for comparing # # two DNA sequences using a ponderation matrix. # # # # Usage : local.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 ### Wi = -1 # penalty for indels ### 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 ### Main ### S = ARGV[0] T = ARGV[1] X = ARGV[2] m = S.length + 1 n = T.length + 1 # Initialisation of ponderation matrix p = PondMatrix.new # Initialisation of array D = Array2D.new for i in 0 ... m D[i,0] = 0 end for j in 0 ... n D[0,j] = 0 end # calculate every other cell for i in 1 ... m for j in 1 ... n D[i,j] = 0 if(D[i,j] < D[i-1,j-1] + p[S[i-1..i-1],T[j-1..j-1]]) D[i,j] = D[i-1,j-1] + p[S[i-1..i-1],T[j-1..j-1]] end if(D[i,j] < D[i,j-1] + Wi) D[i,j] = D[i,j-1] + Wi end if(D[i,j] < D[i-1,j] + Wi) D[i,j] = D[i-1,j] + Wi end end end #Print dynamic table 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("%3d ", D[i,j]) end printf("\n%s ", S[i..i]) end # Print maximal local alignment over S and T. Exercie ; Print all fragments over threshold X without overlaps max = 0 maxI = 0 maxJ =0 for i in 0 ... m for j in 0 ... n if (D[i,j] > max) max = D[i,j] maxI = i maxJ = j end end end if(max == 0) puts("No similar regions between sequences") else i = maxI j = maxJ Sa = String.new("") Ta = String.new("") begin if(D[i,j] == D[i-1,j-1] + p[S[i-1..i-1],T[j-1..j-1]]) Sa << S[i-1..i-1] Ta << T[j-1..j-1] i = i-1 j = j-1 else if(D[i,j] == D[i,j-1] + Wi) Sa << "-" Ta << T[j-1..j-1] j = j-1 else Sa << S[i-1..i-1] Ta << "-" i = i-1 end end end while(D[i,j] > 0) Sa.reverse! Ta.reverse! printf("\nOptimal alignment of score %d:\n", D[maxI,maxJ]) puts(Sa) puts(Ta) puts("\n") end