%IsoFinderSegmentation - DNA segmentation using IsoFinder algorithm*
% 
% Description: IsoFinderSegmentation recursively partitions a sequence using Shannon
% entropy and a threshold calculated based on Monte Carlo simulation. 
% 
% Parameters: IsoFinderSegmentation(dat, target_file, win_size, sizelim_)
%             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.
%             size_lim_ - Optional. Minimum domain size. Should be a multiplication of win_size. Default is: 3008.  
% 
% 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, Domain GC content,
%                     Domain GC content standard deviation
%  
% Examples: IsoFinderSegmentation('Example.txt')
%           IsoFinderSegmentation('Example.txt', 'output.txt')
%           IsoFinderSegmentation('Example.txt', 'output.txt', 32, 3008)
%           IsoFinderSegmentation('Example.txt', 'output.txt', 32, 3008, 0.0002)
%           IsoFinderSegmentation([ 12 13 14 15 18 30 32 31 29 28])
%
% * The script is our interpetation to the algorithm described in:  
%   IsoFinder: computational prediction of isochores in genome sequences. Oliver et al. 2004. Nucleic Acids Res 32:W287-292.
%   Because we could not obtain the source code we had to make several assumption regarding the Monte Carlo simulation. Thereofre, this algorithm
%   is by no mean a representation of the algorithm as it described by the authors (Oliver et al. 2004) and the results may be different from the IsoFinder excutable version.
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Written by Eran Elhaik
% Written in : 07/01/08
% Ver : 1.70
% Website: http://www.igm.jhmi.edu/~eelhaik/Simulation/index.htm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function IsoFinderSegmentation(dat, target_file, win_size, sizelim_)
	
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 = ['IsoFinder_output.txt']; disp(['Output file is ' target_file]); win_size = 32; sizelim_ = 3008; end;
   if (nargin == 2), win_size = 32; sizelim_ = 3008; end;
   if (nargin == 3), sizelim_ = win_size^2*3; 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(['IsoFinderSegmentation Error: Sequence should contain GC counts ranging from 0 to win_size. Exit program.']);
    return;
end;
    
%If sizelim is not a multiplication of win_size adjust it
sizelim = ceil(sizelim_/win_size);
if ~(sizelim==(sizelim_/win_size))
    disp(['IsoFinderSegmentation Warning: It is recommanded to set sizelim as a multiplication of win_size. System would adjust sizelim to ' num2str(sizelim*win_size)]);
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 = CalculateEntropy(dat(fromto(i,1):fromto(i,2))./win_size);

        %If the matrix is not empty
        if (parm(1) ~= 0)
            %If all conditions regarding threshold, min seg, and min sd
            if ( parm(2)>sizelim )
                mat = [mat; fromto(i,1)+parm(1)];
                StopSegmentation = 0;
            end;
        else
            disp(['IsoFinderSegmentation Error: Error on input sequence. ' num2str([i parm]) ]);
            disp(fromto(i,:));
        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

%Calculate statistics for each domain: mean, std
iso_statistics = zeros(2,size(fromto32,1));
for stat = 1:size(fromto32,1)
    seq = dat(fromto(stat,1):fromto(stat,2))/win_size;  %sequence GC%
    iso_statistics(:,stat) = [mean(seq) std(seq)];   %This is std
end;
results = [fromto32 diff(fromto32,1,2)+1 iso_statistics'];

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   CalculateEntropy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out=CalculateEntropy(gc);

%Seq length
len = length(gc);

%Flipping gc
gcr = fliplr(gc);

std_i = 1:len;

%Calculate standard deviation of vector
xbar = cumsum(gc)./std_i;					%The mean GC in every point -> mean(1, 1-2, 1-3, 1-4)
s1 = sqrt((cumsum(gc.^2) - (1:len).*xbar.^2)./(0:(len-1)));
s1(1) = 0;
s1_2 = power(s1,2);

xbarr = fliplr(cumsum(gcr)./(std_i));					%The mean GC in every point
s2 = sqrt(     (cumsum(gcr.^2) - (1:len).*xbarr.^2) ./ (0:(len-1))     );
s2(1) = 0;
s2_2 = fliplr(power(s2,2));

%This is the SD
xstd=((((std_i-1).*s1_2 + (len-std_i).*s2_2)./(len-2)).^0.5) .* ((1./std_i + 1./(len-std_i+1)).^0.5);

%Calculating the t-statistic
t_stat = [];
%t_stat = abs( (xbar-xbarr)./xstd );
t_stat = (xbar-xbarr)./xstd;

%finding t max
[t_max_value t_max_loc] = max(abs(t_stat));
[t_max_value t_max_loc] = max(t_stat);

%Building an index of 3k numbers from 1 to size of gc. 
L0 = 94;

%If t_max_loc is within boubderies of L0 size.
%If it is not - it will not be possible to calculate the GC in L0 size
%window
if ((t_max_loc > L0) & (len-t_max_loc> L0))

    %Calcualte indexes for L0bp windows from right and left of t_max
    L0left = gc(1:L0:t_max_loc);    %Calculate the values only in L0+1 positions
    L0right = gc(t_max_loc:L0:len);

    %Creating temp vectors for left and right parts of t max
    gc_left = gc(1:L0*(length(L0left)-1));
    gc_right = gc(t_max_loc:t_max_loc-1+(L0*(length(L0right)-1)));

    %Calculating the mean for every L0 window     
    mean_gc_left = mean(reshape(gc_left,length(L0left)-1, L0)');
    mean_gc_right = mean(reshape(gc_right,length(L0right)-1, L0)');

    %Finally, we check whether this significance
    %exceeds a selected threshold P0, usually taken to be 95%.
    
    %Make sure there are enough cells to use ttest
    if ((length(mean_gc_left)==1) & (length(mean_gc_right)==1))
        out=[0 0]; 
    elseif ((sum(mean_gc_left)==0) & (sum(mean_gc_right)==0))
        out=[0 0]; 
    else
    
        %Student ttest between right and left means    
        [h p ci stats]=ttest2(mean_gc_left, mean_gc_right);

	%Start Monte Carlo simulations
	
	%Calculate the beta incomplete function
	N=len;
	v=N-2;
	Ix = v/(v+power(stats.tstat,2));
	Iz = 0.4*v;
	Iw = 0.4;
	lambda = (4.19*log(N))-11.54;
	
	I = betainc(Ix, Iz, Iw);
	P_t = (1-I)^lambda;
	
	if P_t>0.95
           out=[1 t_max_loc];
        else
           out=[2 0]; 
        end;
    end;
else
    %If i will send t_max_loc here it will result in endless loop
    out=[3 0]; 
end;

