%CostantiniSegmentation - DNA segmentation using sliding window as proposed by Costantini et al. (2006)*
% 
% Description: CostantiniSegmentation uses a sliding window of size 100,000bp to partition a genomic sequence. 
% 
% Parameters: CostantiniSegmentation(dat, target_file, win_size)
%             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 lengt.
%  
% Examples: CostantiniSegmentation('Example.txt')
%           CostantiniSegmentation('Example.txt', 'output.txt')
%           CostantiniSegmentation('Example.txt', 'output.txt', 32)
%           CostantiniSegmentation([ 12 13 14 15 18 30 32 31 29 28])
%
% * Because the algorithm code is unavailable we coded the algorithm based on the description in 
%   An isochore map of human chromosomes. Costantini et al. 2006. Genome Res. 16:536-541.
%   We note that this algorithm represents our understanding to the method proposed by Costantini et al. 2006
%   and results may be different from the results published by the authors.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Written by Eran Elhaik
% Written in : 07/01/08
% Ver : 1.00
% Website: http://www.igm.jhmi.edu/~eelhaik/Simulation/index.htm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CostantiniSegmentation(dat, target_file, win_size)

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

threshold = 0.01; %value set according to paper

%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(['CostantiniSegmentation Error: Sequence should contain GC counts ranging from 0 to win_size. Exit program.']);
    return;
end;
    
window_size = floor(100000/win_size);

borders = 1:window_size:length(dat);

len = length(borders);

b = reshape(dat(1:borders(len)-1),window_size,len-1);

c = std(b);

%Entering a first number to avoid a frame shift caused by diff
d = [9 abs(diff(c))];

% Check treshold 
segments = find (d > threshold);

if length(segments) > 1
   segments = segments * window_size;
   segments = [1 segments*win_size length(dat)*win_size];
else
   segments = [1 length(dat)*win_size];
end;

segments=segments';

%Finalize vector
fromto = [segments(1:1:end-1) segments(2:1:end) segments(2:1:end)-segments(1:1:end-1)]; %Convert to 3 columns format

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

