Code covered by the BSD License  

Highlights from
Calculate number of cycles in genetic 2-Break distance problem

image thumbnail

Calculate number of cycles in genetic 2-Break distance problem

by

 

25 Jan 2013 (Updated )

Given two genomes, find the shortest sequence of 2-Breaks transforming one genome into another

Dist_2break()
function Dist_2break()
  P = [1, 2, 3, 4];
  Q = [1, -3, -2, 4];  
  
  
  len = numel(P);
  for i = 1:len
      Qp(i) = find(abs(P) == abs(Q(i)));
  end
  Qp = Qp .* sign(P(Qp)) .* sign(Q);
  %% update map table for P
  % - for tail, + for head
  for i = 1:len
      tP(i, 1) = sign(P(i))*i;
      tP(i, 2) = -sign(P(mod(i, len)+1))*(mod(i, len)+1);
  end
  %% update map table for Q
  for i = 1:len
      tQ(i, 1) = sign(P(abs(Qp(i)))) * Qp(i);
      tQ(i, 2) = -sign(P(abs(Qp(mod(i,len)+1)))) * Qp(mod(i,len)+1);
  end
  
  %% get cycles
  cycle_count = 0;  
  i_p = 1;
  j_p = 1;
  c_start = [];
  for i = 1:len
      [c_start, c_end, tP, tQ] = updateCycle(tP, tQ, i_p, j_p, c_start);         
      if(c_start == c_end)
          cycle_count = cycle_count + 1;
          c_start = []; 
          i_p = 1;  j_p = 1;
      else     
      [i_p, j_p] = findNeighborIndex(tP, c_end); 
      end
  end
  cycle_count
end
 
function [c_start, c_end, tP, tQ] = updateCycle(tP, tQ, i_p, j_p, c_start)      
      [i_ref, j_ref] = findIndex(tQ, tP(i_p, j_p));
      %% update current cycle end
      if(isempty(c_start))
        c_start = tP(i_p, 3 - j_p); 
      end
      c_end = tQ(i_ref, 3 - j_ref); 
      %% remove processed row.
      tQ(i_ref, :) = [];
      tP(i_p, :) = [];
end

% find 2D index in tQ for the given Value
function [i, j] = findIndex(A, val)
    [i, j] = ind2sub(size(A), find(A(:) == val));
end

% find 2D index in tQ for the given Value
function [i, j] = findNeighborIndex(A, val)
    [i, j] = findIndex(A, val);
    j = 3 - j;
end
  
  

Contact us