disp('Running MATLAB script file LongDat1fig6.m') ;
%
%    Makes Figure 6 for the Longitudinal Data Analysis paper
% 
% angle boundaries are set by examing the Sizer plot only
% G1 = [2.04 3.74); S = [3.74, 4.48); G2 = [4.48 5.72); M = [5.72 6.280) and [0, 0.83); MG1 = [0.83 2.04);

% 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 200 genes
distance=sqrt(proj_sin.^2 + proj_cos.^2);
[dist_order, vind]=sort(distance);
vind_200=vind(4290:4489); % get the index for top 200 genes
vind_no=vind(1:4289)'; % get the index for the rest of genes.
angle_200=angle(vind_200); % 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_200)
    if (angle_200(i) >= G1_low) & (angle_200(i) < S_low)
        %angle_G1 = [angle_G1;angle_regu(i)];
        vind_G1=[vind_G1;vind_200(i)];
    elseif (angle_200(i) >= S_low) & (angle_200(i) < G2_low)
        %angle_S = [angle_S;angle_regu(i)];
        vind_S=[vind_S;vind_200(i)];
    elseif (angle_200(i) >= G2_low) & (angle_200(i) < M_low)
        %angle_G2 = [angle_G2;angle_regu(i)];
        vind_G2=[vind_G2;vind_200(i)];
    elseif ((angle_200(i) >= M_low) & (angle_200(i) <= 2*pi)) | ((angle_200(i) >= 0) & (angle_200(i) < MG1_low))
        %angle_M = [angle_M;angle_regu(i)];
        vind_M=[vind_M;vind_200(i)];  
    elseif (angle_200(i) >= MG1_low) & (angle_200(i) < G1_low)
        %angle_MG1 = [angle_MG1;angle_regu(i)];
        vind_MG1=[vind_MG1;vind_200(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));

% 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];

proj_sin_order=proj_sin(vind_order);
proj_cos_order=proj_cos(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 ;

plot(proj_cos_order(1), proj_sin_order(1),'+','Color', color_order(1,:));
hold on;
for idat=2:length(vind_order)
    plot(proj_cos_order(idat), proj_sin_order(idat),'+','Color', color_order(idat,:));
end
plot([0 2*cos(G1_low)], [0 2*sin(G1_low)], 'm-');
plot([0 2*cos(S_low)], [0 2*sin(S_low)], 'm-');
plot([0 3*cos(G2_low)], [0 3*sin(G2_low)], 'm-');
plot([0 2*cos(M_low)], [0 2*sin(M_low)], 'm-');
plot([0 2*cos(MG1_low)], [0 2*sin(MG1_low)], 'm-');
hold off;
%title('scatterplot of ProjSin vs ProjCos for all genes (4489)');
%xlabel('ProjCos');
%ylabel('ProjSin');
axis([-2 2 -2 2]);
set(gca,'XTickLabelMode','manual') ;
set(gca,'YTickLabelMode','manual') ;




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