filename    = 'Katrina_Latest_Statistical_Analysis_Normalized_Data_20111101.xls';
sheetname   = 'Master_Sheet';
strain1     = 'WildType';
strain2     = 'dCIN5';
strain3     = 'dGLN3';
strain4     = 'dHMO1';
strain5     = 'dZAP1';

[log2FC,b]=xlsread(filename,sheetname);
a = log2FC(:,2:end);
%
% a = numerical data from sheet in matrix form
% b = character data from sheet in cell array form
%
%
% Set up indices to access the data.  These are the columns
%   of the a matrix corresponding to each time point.
%   This is not automated and must be edited by the user
%

ind15  = [1,2,3,4];
ind30  = [5,6,7,8,9];
ind60  = [10,11,12,13];
ind90  = [14,15,16,17,18];
ind120 = [19,20,21,22,23];

ind152  = [33,34,35,36];
ind302  = [37,38,39,40];
ind602  = [41,42,43,44];
ind902  = [45,46,47,48];
ind1202 = [49,50,51,52];

ind153  = [62,63,64,65];
ind303  = [66,67,68,69];
ind603  = [70,71,72,73];
ind903  = [74,75,76,77];
ind1203 = [78,79,80,81];

ind154  = [91,92,93,94];
ind304  = [95,96,97,98];
ind604  = [99,100,101,102];
ind904  = [103,104,105,106];
ind1204 = [107,108,109,110];

ind155  = [120,121,122,123];
ind305  = [124,125,126,127];
ind605  = [128,129,130,131];
ind905  = [132,133,134,135];
ind1205 = [136,137,138,139];

p = 25;
q = 20;

alpha = 0.05;  % significance level for the tests

n152 = length(ind152);
n302 = length(ind302);
n602 = length(ind602);
n902 = length(ind902);
n1202 = length(ind1202);

n153 = length(ind153);
n303 = length(ind303);
n603 = length(ind603);
n903 = length(ind903);
n1203 = length(ind1203);

n154 = length(ind154);
n304 = length(ind304);
n604 = length(ind604);
n904 = length(ind904);
n1204 = length(ind1204);

n155 = length(ind155);
n305 = length(ind305);
n605 = length(ind605);
n905 = length(ind905);
n1205 = length(ind1205);

indx2   = [ind152,ind302,ind602,ind902,ind1202];
indx3   = [ind153,ind303,ind603,ind903,ind1203];
indx4   = [ind154,ind304,ind604,ind904,ind1204];
indx5   = [ind155,ind305,ind605,ind905,ind1205];

n = length(a(:,1));

%
% out_data includes
%
out_data = zeros(n,12);
nsig = 0;

% % These next three lines of code enable a movie to be written of the plots
% % that will be produce.
% figure(1)
% set(gcf, 'DoubleBuffer','on');
% yourMovieFileName = ['all_strain_comparison_p_value_less_than_0.05_Cinepak.avi'];

for ii=1:n
    I = find(~isnan(a(ii,1:23)));
    ind15 = I(find(I>=1&I<=4));
    ind30 = I(find(I>=5&I<=9));
    ind60 = I(find(I>=10&I<=13));
    ind90 = I(find(I>=14&I<=18));
    ind120 = I(find(I>=19&I<=23));

    indx    = [ind15,ind30,ind60,ind90,ind120];
    indX    = find(indx);

    if length(indx)<23
        ind15 = find(ind15);
        ind30 = ind30 - ind15(end);
        ind60 = ind60 - ind30(1) - 1;
        ind90 = ind90 - ind60(1);
        ind120 = ind120-ind90(1);
    end

    n15 = length(ind15);
    n30 = length(ind30);
    n60 = length(ind60);
    n90 = length(ind90);
    n120 = length(ind120);

    i1 = indX(end);
    i2 = indx2(1);

    ind152m = ind152-i2+1+i1;
    ind302m = ind302-i2+1+i1;
    ind602m = ind602-i2+1+i1;
    ind902m = ind902-i2+1+i1;
    ind1202m = ind1202-i2+1+i1;

    i3 = ind1202m(end);
    i4 = indx3(1);

    ind153m = ind153-i4+1+i3;
    ind303m = ind303-i4+1+i3;
    ind603m = ind603-i4+1+i3;
    ind903m = ind903-i4+1+i3;
    ind1203m = ind1203-i4+1+i3;

    i5 = ind1203m(end);
    i6 = indx4(1);

    ind154m = ind154-i6+1+i5;
    ind304m = ind304-i6+1+i5;
    ind604m = ind604-i6+1+i5;
    ind904m = ind904-i6+1+i5;
    ind1204m = ind1204-i6+1+i5;

    i7 = ind1204m(end);
    i8 = indx5(1);

    ind155m = ind155-i8+1+i7;
    ind305m = ind305-i8+1+i7;
    ind605m = ind605-i8+1+i7;
    ind905m = ind905-i8+1+i7;
    ind1205m = ind1205-i8+1+i7;

    t = [ones(n15,1)*15;ones(n30,1)*30;ones(n60,1)*60;ones(n90,1)*90;ones(n120,1)*120;...
        ones(n152,1)*15;ones(n302,1)*30;ones(n602,1)*60;ones(n902,1)*90;ones(n1202,1)*120;...
        ones(n153,1)*15;ones(n303,1)*30;ones(n603,1)*60;ones(n903,1)*90;ones(n1203,1)*120;...
        ones(n154,1)*15;ones(n304,1)*30;ones(n604,1)*60;ones(n904,1)*90;ones(n1204,1)*120;...
        ones(n155,1)*15;ones(n305,1)*30;ones(n605,1)*60;ones(n905,1)*90;ones(n1205,1)*120];

    t2 = [min(t):5:max(t)]';

    N = length(indx)+length(indx2)+length(indx3)+length(indx4)+length(indx5);

    X   = zeros(N,p);
    Xh  = zeros(N,p-q);

    Xh(ind15,1) = 1;
    Xh(ind30,2) = 1;
    Xh(ind60,3) = 1;
    Xh(ind90,4) = 1;
    Xh(ind120,5) = 1;
    Xh(ind152m,1) = 1;
    Xh(ind302m,2) = 1;
    Xh(ind602m,3) = 1;
    Xh(ind902m,4) = 1;
    Xh(ind1202m,5) = 1;
    Xh(ind153m,1) = 1;
    Xh(ind303m,2) = 1;
    Xh(ind603m,3) = 1;
    Xh(ind903m,4) = 1;
    Xh(ind1203m,5) = 1;
    Xh(ind154m,1) = 1;
    Xh(ind304m,2) = 1;
    Xh(ind604m,3) = 1;
    Xh(ind904m,4) = 1;
    Xh(ind1204m,5) = 1;
    Xh(ind155m,1) = 1;
    Xh(ind305m,2) = 1;
    Xh(ind605m,3) = 1;
    Xh(ind905m,4) = 1;
    Xh(ind1205m,5) = 1;

    X(ind15,1) = 1;
    X(ind30,2) = 1;
    X(ind60,3) = 1;
    X(ind90,4) = 1;
    X(ind120,5) = 1;
    X(ind152m,6) = 1;
    X(ind302m,7) = 1;
    X(ind602m,8) = 1;
    X(ind902m,9) = 1;
    X(ind1202m,10) = 1;
    X(ind153m,11) = 1;
    X(ind303m,12) = 1;
    X(ind603m,13) = 1;
    X(ind903m,14) = 1;
    X(ind1203m,15) = 1;
    X(ind154m,16) = 1;
    X(ind304m,17) = 1;
    X(ind604m,18) = 1;
    X(ind904m,19) = 1;
    X(ind1204m,20) = 1;
    X(ind155m,21) = 1;
    X(ind305m,22) = 1;
    X(ind605m,23) = 1;
    X(ind905m,24) = 1;
    X(ind1205m,25) = 1;


    Y = [a(ii,indx)';a(ii,indx2)';a(ii,indx3)';a(ii,indx4)';a(ii,indx5)'];

    beta    = X\Y;
    betah   = Xh\Y;

    out_data(ii,1:25) = beta';
    out_data(ii,26:30) = betah';

    s2 = (1/N)*(Y-X*beta)'*(Y-X*beta);
    s2h = (1/N)*(Y-Xh*betah)'*(Y-Xh*betah);

    W       = ((N-p)/q)*((s2h-s2)/s2);
    pval    = 1-fcdf(W,q,N-p);

    out_data(ii,31) = W;
    out_data(ii,32) = pval;
    out_data(ii,33) = s2*N;
    out_data(ii,34) = s2h*N;

    Ym  = X*beta;
    Yh  = Xh*betah;

%         Turns on the capability to pause at a graph.
    pause on
    if pval<0.05
        figure(1),plot(t,Y,'bo'),hold on,plot(t,Ym,'r+','LineWidth',2),plot(t,Yh,'k^','LineWidth',2),...
            title([b{ii+1,1} '/' b{ii+1,2}],'FontSize',12),...
            xlabel({'time, minutes', '' , ['   W = ',num2str(W),' pval = ',num2str(pval)]}),...
            ylabel('log fold change'),legend('data','model'),drawnow
        hold off
        nsig = nsig + 1;
        [ii,nsig]
        drawnow
% %         Save as a .jpg
%         saveas(figure(1),[b{ii+1,1} '.jpg'])
        % %         Save as a .fig
        %         saveas(figure(1),[b{ii+1,1} '.fig'])
        %         Allows you to pause at a graph for as long as you do not press a key on the keyboard.
        pause
        % %         Captures the figure being currently displayed to put into the movie.
        %         movieArray(nsig) = getframe(gcf);
    end

end

% %THis complies your movie.
% movie2avi(movieArray,yourMovieFileName,'compression','Cinepak','fps',5,'quality',100);

p1 = out_data(:,32);
q1 = (alpha/n)*[1:n]';

[ps,is] = sort(p1);

tests   = ps<=q1;
itests  = find(tests==1);
itestm  = max(itests);

qs = zeros(size(q1));
rs = qs;
for ii = 1:n
    qs(is(ii)) = ps(ii)*n/ii;
    rs(is(ii)) = ii;
end

signif                  = zeros(n,1);
signif(is(1:itestm))    = 1;

out_data(:,35) = qs;
out_data(:,36) = signif;
out_data(:,37) = rs;

eval(['save ' strain1 '_' strain2 '_' strain3 '_' strain4 '_' strain5 '_out_data out_data;']);

out_data_cells{1,1} = 'Systematic Name';
out_data_cells{1,2} = 'Standard Name';
out_data_cells{1,3} = 't15';
out_data_cells{1,4} = 't30';
out_data_cells{1,5} = 't60';
out_data_cells{1,6} = 't90';
out_data_cells{1,7} = 't120';
out_data_cells{1,8} = 't15';
out_data_cells{1,9} = 't30';
out_data_cells{1,10} = 't60';
out_data_cells{1,11} = 't90';
out_data_cells{1,12} = 't120';
out_data_cells{1,13} = 't15';
out_data_cells{1,14} = 't30';
out_data_cells{1,15} = 't60';
out_data_cells{1,16} = 't90';
out_data_cells{1,17} = 't120';
out_data_cells{1,18} = 't15';
out_data_cells{1,19} = 't30';
out_data_cells{1,20} = 't60';
out_data_cells{1,21} = 't90';
out_data_cells{1,22} = 't120';
out_data_cells{1,23} = 't15';
out_data_cells{1,24} = 't30';
out_data_cells{1,25} = 't60';
out_data_cells{1,26} = 't90';
out_data_cells{1,27} = 't120';
out_data_cells{1,28} = 't15';
out_data_cells{1,29} = 't30';
out_data_cells{1,30} = 't60';
out_data_cells{1,31} = 't90';
out_data_cells{1,32} = 't120';
out_data_cells{1,33} = 'f stat';
out_data_cells{1,34} = 'p val';
out_data_cells{1,35} = 'SS full';
out_data_cells{1,36} = 'SS hyp';
out_data_cells{1,37} = 'B&H comps';
out_data_cells{1,38} = '? signif ?';
out_data_cells{1,39} = 'B&H rank';
for ii = 1:n
    out_data_cells{1+ii,1} = b{ii+1,1};
    out_data_cells{1+ii,2} = b{ii+1,2};
    for jj = 1:37
        out_data_cells{1+ii,2+jj} = out_data(ii,jj);
    end
end
xlswrite([strain1 '_' strain2 '_' strain3 '_' strain4 '_' strain5 '_out_data.xls'],out_data_cells)