Code covered by the BSD License

# Calculate number of cycles in genetic 2-Break distance problem

### Shuang Wang (view profile)

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

```