disp('Running MATLAB script file LongDat1fig5.m') ;
%
%    Makes Figure 5 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
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;

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

%paramstruct_sizer=struct('ibdryadj', 2, 'iout', 2);
paramstruct_sizer=struct('ibdryadj', 2, 'iout', 2, 'fhmin',10^(-1.5), ...
                                                   'fhmax',10^(0.2), ...
                                                   'shmin',10^(-1.5), ...
                                                   'shmax',10^(0.2));
makeplot=sizerSM(angle_200, paramstruct_sizer);
subplot(3,1,1),
vax=axis;
hold on;
plot([G1_low G1_low],[vax(3) vax(4)],'k-');
plot([S_low S_low],[vax(3) vax(4)],'k-');
plot([G2_low G2_low],[vax(3) vax(4)],'k-');
plot([M_low M_low],[vax(3) vax(4)],'k-');
plot([MG1_low MG1_low],[vax(3) vax(4)],'k-');
hold off;
subplot(3,1,2),
vax=axis;
hold on;
plot([G1_low G1_low],[vax(3) vax(4)],'k-');
plot([S_low S_low],[vax(3) vax(4)],'k-');
plot([G2_low G2_low],[vax(3) vax(4)],'k-');
plot([M_low M_low],[vax(3) vax(4)],'k-');
plot([MG1_low MG1_low],[vax(3) vax(4)],'k-');
hold off;
subplot(3,1,3),
vax=axis;
hold on;
plot([G1_low G1_low],[vax(3) vax(4)],'k-');
plot([S_low S_low],[vax(3) vax(4)],'k-');
plot([G2_low G2_low],[vax(3) vax(4)],'k-');
plot([M_low M_low],[vax(3) vax(4)],'k-');
plot([MG1_low MG1_low],[vax(3) vax(4)],'k-');
hold off;

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