disp('Running MATLAB script file LongDat1fig8.m') ;
%
%    Makes Figure 8 for the Longitudinal Data Analysis paper
% 

% choose threshold value for identifying periodic genes
thre=612;

% Input data
fid=fopen('spellman_alpha_complete.txt', 'rt');
genevec=fscanf(fid, '%f');
mdata=reshape(genevec, [18, 4489]);  % dimention = # of time points

% center the data, and compute the mean residuals.
vcenter=mean(mdata,2);
mresid=mdata-vec2matSM(vcenter,4489);

% get projections of the mean residuals onto cos(4pi) and sin(4pi) spaces.
d=18;
j=1:d;
w=4*pi/d;
cosX=[cos(w*j)]';
sinX=[sin(w*j)]';
proj_cos=inv(cosX'*cosX)*cosX'*mresid;
proj_sin=inv(sinX'*sinX)*sinX'*mresid;

% compute the angles ([-pi, pi])of each gene on the scatterplot of the two projections.
angle=atan2(proj_sin, proj_cos)'; %angle is 4489*1 column vector

% change angles in (-pi, 0) to (pi, 2*pi), so angles on [0 2*pi]
for i=1:length(angle)
    if angle(i)<0
        angle(i)=2*pi+angle(i);
    end
end

% sort the distances to the origin of each point in the scatterplot 
% order them in ascending order
% get the index for the top thre genes
distance=sqrt(proj_sin.^2 + proj_cos.^2);
[dist_order, vind]=sort(distance);

index_thre=4489-thre+1;
vind_thre=vind(index_thre:4489); % get the index for top thre genes
vind_no=vind(1:index_thre-1)'; % get the index for the rest of genes.
angle_thre=angle(vind_thre); % get the top # of genes which we think are cell cycle-regulated

% assign the angle boundaries for five phases
G1_low = 2.04;
S_low = 3.74;
G2_low = 4.48;
M_low = 5.72;
MG1_low = 0.83;

% classify top 200 genes into five phase groups according to paramstruct.
vind_G1=0;
vind_S=0;
vind_G2=0;
vind_M=0;
vind_MG1=0;
for i=1:length(angle_thre)
    if (angle_thre(i) >= G1_low) & (angle_thre(i) < S_low)
        vind_G1=[vind_G1;vind_thre(i)];
    elseif (angle_thre(i) >= S_low) & (angle_thre(i) < G2_low)
        vind_S=[vind_S;vind_thre(i)];
    elseif (angle_thre(i) >= G2_low) & (angle_thre(i) < M_low)
        vind_G2=[vind_G2;vind_thre(i)];
    elseif ((angle_thre(i) >= M_low) & (angle_thre(i) <= 2*pi)) | ((angle_thre(i) >= 0) & (angle_thre(i) < MG1_low))
        vind_M=[vind_M;vind_thre(i)];  
    elseif (angle_thre(i) >= MG1_low) & (angle_thre(i) < G1_low)
        vind_MG1=[vind_MG1;vind_thre(i)];
    end
end

vind_G1=vind_G1(2:length(vind_G1));
vind_S=vind_S(2:length(vind_S));
vind_G2=vind_G2(2:length(vind_G2));
vind_M=vind_M(2:length(vind_M));
vind_MG1=vind_MG1(2:length(vind_MG1));

count_G1=length(vind_G1);
count_S=length(vind_S);
count_G2=length(vind_G2);
count_M=length(vind_M);
count_MG1=length(vind_MG1);


% get ordered proj_cos and proj_sin for all genes.
% order is nonperiodic->G1->S->G2->M->M/G1
vind_order=[vind_no;vind_G1;vind_S;vind_G2;vind_M;vind_MG1];

mdata_order=mdata(:,vind_order);

for i=1:length(vind_no)
icolor_no(i,:)=[0 0 0];  % angle_no is black
end
for i=1:length(vind_G1)
icolor_G1(i,:)=[1 0 0];  % G1 phase is red
end
for i=1:length(vind_S)
icolor_S(i,:)=[0 1 0];  % S phase is green
end
for i=1:length(vind_G2)
icolor_G2(i,:)=[0 0 1];  % G2 phase is blue
end
for i=1:length(vind_M)
icolor_M(i,:)=[1 1 0];  % M phase is yellow
end
for i=1:length(vind_MG1)
icolor_MG1(i,:)=[0 1 1];  % M/G1 phase is cyan
end

%get color map for angle_order
color_order = [icolor_no;icolor_G1;icolor_S;icolor_G2;icolor_M;icolor_MG1];

%  make main graphic
%
figure(1) ;
clf ;

% top thre genes
index_thre=4489-thre+1;
ah1 = subplot(3,2,1);
    paramstruct_curv_200=struct('viout', 1, 'vipcplot', 0, 'vicolplot', 1, 'icolor', color_order(index_thre:4489,:), 'isingleaxis', 1, 'iscreenwrite', 1);
    CurvDatSM(mdata_order(:,index_thre:4489),paramstruct_curv_200); 
    axis([1 18 -3 5]);
    title('') ;
    set(ah1,'XTickLabelMode','manual') ;
    set(ah1,'XTickLabel',[]) ;
    set(ah1,'YTickLabelMode','manual') ;
    set(ah1,'yTickLabel',[]) ;
    vchil = get(ah1,'Children') ;

% G1 group
index_G1=index_thre + count_G1 - 1;
ah2 = subplot(3,2,2);
    paramstruct_curv_G1=struct('viout', 1, 'vipcplot', [0], 'vicolplot', [1], 'icolor', color_order(index_thre:index_G1,:), 'isingleaxis', 1, 'iscreenwrite', 1);
    CurvDatSM(mdata_order(:,index_thre:index_G1),paramstruct_curv_G1); 
    axis([1 18 -3 5]);
    title('') ;
    set(ah2,'XTickLabelMode','manual') ;
    set(ah2,'XTickLabel',[]) ;
    set(ah2,'YTickLabelMode','manual') ;
    set(ah2,'yTickLabel',[]) ;
    vchil = get(ah2,'Children') ;
    %set(vchil(1),'String',[]) ;
    
% S group
index_S=index_G1 + count_S;
ah3 = subplot(3,2,3);
    paramstruct_curv_S=struct('viout', 1, 'vipcplot', [0], 'vicolplot', [1], 'icolor', color_order((index_G1+1):index_S,:), 'isingleaxis', 1, 'iscreenwrite', 1);
    CurvDatSM(mdata_order(:,(index_G1+1):index_S),paramstruct_curv_S); 
    axis([1 18 -3 5]);
    title('') ;
    set(ah3,'XTickLabelMode','manual') ;
    set(ah3,'XTickLabel',[]) ;
    set(ah3,'YTickLabelMode','manual') ;
    set(ah3,'yTickLabel',[]) ;
    vchil = get(ah3,'Children') ;
    %set(vchil(1),'String',[]) ;

% G2 group
index_G2=index_S + count_G2;
ah4 = subplot(3,2,4);
    paramstruct_curv_G2=struct('viout', 1, 'vipcplot', [0], 'vicolplot', [1], 'icolor', color_order((index_S+1):index_G2,:), 'isingleaxis', 1, 'iscreenwrite', 1);
    CurvDatSM(mdata_order(:,(index_S+1):index_G2),paramstruct_curv_G2); 
    axis([1 18 -3 5]);
    title('') ;
    set(ah4,'XTickLabelMode','manual') ;
    set(ah4,'XTickLabel',[]) ;
    set(ah4,'YTickLabelMode','manual') ;
    set(ah4,'yTickLabel',[]) ;
    vchil = get(ah4,'Children') ;
    %set(vchil(1),'String',[]) ;

% M group
index_M=index_G2 + count_M;
ah5 = subplot(3,2,5);
    paramstruct_curv_M=struct('viout', 1, 'vipcplot', [0], 'vicolplot', [1], 'icolor', color_order((index_G2+1):index_M,:), 'isingleaxis', 1, 'iscreenwrite', 1);;
    CurvDatSM(mdata_order(:,(index_G2+1):index_M),paramstruct_curv_M); 
    axis([1 18 -3 5]);
    title('') ;
  set(ah5,'XTickLabelMode','manual') ;
  set(ah5,'XTickLabel',[]) ;
  set(ah5,'YTickLabelMode','manual') ;
  set(ah5,'yTickLabel',[]) ;
  vchil = get(ah5,'Children') ;
  %set(vchil(1),'String',[]) ;

% MG1 group
index_MG1=index_M + count_MG1;
ah6 = subplot(3,2,6);
    paramstruct_curv_MG1=struct('viout', 1, 'vipcplot', [0], 'vicolplot', [1], 'icolor', color_order((index_M+1):index_MG1,:), 'isingleaxis', 1, 'iscreenwrite', 1);
    CurvDatSM(mdata_order(:,(index_M+1):index_MG1),paramstruct_curv_MG1); 
    axis([1 18 -3 5]);
    title('') ;
  set(ah6,'XTickLabelMode','manual') ;
  set(ah6,'XTickLabel',[]) ;
  set(ah6,'YTickLabelMode','manual') ;
  set(ah6,'yTickLabel',[]) ;
  vchil = get(ah6,'Children') ;
  %set(vchil(1),'String',[]) ;


 pstr = 'LongDat1fig9a' ;
 orient portrait ;
 set(gcf,'PaperSize',[6.0, 5.0]) ; 
 set(gcf,'PaperPosition',[0.25, 0.25, 5.0, 4.5]) ; 
print('-dpsc', [pstr '.ps']) ;