% ChiSegmentation - DNA segmentation using CHI 
% 
% Description: ChiSegmentation recursively partitions a sequence using CHI method, proposed by Nekrutenko and Li (2000)
% 
% Parameters: ChiSegmentation(dat, target_file, win_size, sizelim_, threshold)
%             dat - A txt file that contains the number of GC nucleotides
%                   per a certain window size (default is 32bp). The script accepts either the filename or a vector. 
%		    To create a file with the format, use the bash script: CalculateGCWindow.n.	
%             target_file -  Optional. Output filename. 
%             win_size - Optional. The window for which nucleotides were
%                        counted. Must be the same value as used to create the file. Default is 32bp.
% 
% Input: The program accepts both filename or a vector of GC content counts
%        per window size.
%           Example of input file of 32bp: 16 17 1 2 3 32 14 5 6
%           Example of input file of 1bp:  1  0  0 1 1 1  0  0 0 
% 
% Output file format: From base, To base, Domain length
%  
% Examples: ChiSegmentation('Example.txt')
%           ChiSegmentation('Example.txt', 'output.txt')
%           ChiSegmentation('Example.txt', 'output.txt', 32)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Written by Eran Elhaik
% Written in : 07/01/08
% Ver : 1.20
% Website: http://www.igm.jhmi.edu/~eelhaik/Simulation/index.htm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ChiSegmentation(dat, target_file, win_size)

toplim = 0.05;
sizelim=100;

if (nargin>0)
   if (nargin==0), disp('Not enough parameters! Format is: Source filename, Target filename (optional), window size (optional), minimum domain size (optional)'); return; end;
   if (nargin == 1), target_file = ['CHI_output.txt']; disp(['Output file is ' target_file]); win_size = 32; end;
   if (nargin == 2), win_size = 32; end;
end

%If the dat file is a reference to a file
if ~isnumeric(dat)
    dat=load(dat);  %If this is a file link - load the file
end; 

%Create a row vector of data
dat=reshape(dat, [1 length(dat) ]);  

%Quality control of input sequence
if length([find(dat<0) find(dat>win_size)])
    disp(['ChiSegmentation Error: Sequence should contain GC counts ranging from 0 to win_size. Exit program.']);
    return;
end;

StopSegmentation = 0;
mat=[1; length(dat)];

%Start recursive segmentation of dat
while StopSegmentation<1,
    StopSegmentation=1;

    %Get fromto sizes (2 vecors with from and to)
    fromto = [mat(1:1:end-1) mat(2:1:end)];

    %Run on the length of the vecotr and for each border
    for i=1:size(fromto,1)
        parm = CalculateCHI(dat(fromto(i,1):fromto(i,2)), win_size);


        %If the matrix is not empty
        if ((size(parm,2)) ~= 0)
            if (parm(2)>toplim & parm(1)>sizelim & fromto(i,2)-parm(1)>sizelim & fromto(i,2)-(fromto(i,1)+parm(1)) > sizelim)
               mat=[mat; fromto(i,1)+parm(1)];
               stop=0;
            end;
        end;

    end;
    mat=sort(mat);
end;


%Finalize vector
fromto = [mat(1:1:end-1) mat(2:1:end)];         %Convert to 2 columns format
fromto32 = fromto*win_size; fromto32(:,1)=fromto32(:,1)+1; fromto32(1) = 1;    %Convert to Xbp format

results = [fromto32 diff(fromto32)'];

%Write results to file
fid = fopen(target_file, 'w+');
if fid
    fprintf(fid, '%d %d %d \n', results');
    fclose(fid);
else
    disp(['ChiSegmentation Error: Cannnot open output file ' fout]);
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   Calculate CHI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out=CalculateCHI(gc, win_size);

M1 = gc; %dat(fromto(i,1):fromto(i,2));

%Window size in which we calculate the index
N = 2;  

M1_index=1:N:length(M1);

vec = [];

%Run on all groups and calculate index
for j=1:length(M1_index)-1

    M2 = M1((M1_index(j)):(M1_index(j+1)));
        
    %Check that this is necessary and not affects the final result
    GC = M2/win_size;

    %2 unnecessary parameters
    n = length(M2);
    l = 1;

    %Looking in a certain window size, we don't use windows
    P = mean(GC)/l;

    %Calculate dinominator
    SE = power(((P*(1-P))/l),0.5);

    i_index=(2:n);
    j_index=(1:n-1);

    delta_gc1 = (1/(n-1));
    delta_gc2 = sum(abs(GC(i_index)-GC(j_index)));
    delta_gc = delta_gc1 * delta_gc2;

    if (SE == 0)
        Hgc = 0;
    else
        Hgc = delta_gc / SE;
    end;
    
    vec = [vec repmat(Hgc,1,N)];
end;


[Vmax Dmax] = max(vec(1:end-1));

out=[Dmax Vmax];
