disp('Running MATLAB script file LongDat1fig7.m') ;
%
%    Makes Figure 7 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 number of genes in each phase group.
count_all=length(vind_200);
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);

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

% Kde plot of top # angles for specified global bandwidth: h.
h = 0.266;
paramstruct_kde=struct('vh', h, 'vxgrid', [0;2*pi], 'ibdryadj', 2);
[kde_all, xgrid_all, mker_all]=kdeSM(angle_200, paramstruct_kde);
[kde_G1, xgrid_G1, mker_G1]=kdeSM(angle(vind_G1), paramstruct_kde);
[kde_S, xgrid_S, mker_S]=kdeSM(angle(vind_S), paramstruct_kde);
[kde_G2, xgrid_G2, mker_G2]=kdeSM(angle(vind_G2), paramstruct_kde);
[kde_M, xgrid_M, mker_M]=kdeSM(angle(vind_M), paramstruct_kde);
[kde_MG1, xgrid_MG1, mker_MG1]=kdeSM(angle(vind_MG1), paramstruct_kde);

plot(xgrid_all, kde_all, 'k--'); % plot kde for all top # gene angles in black
hold on;
plot(xgrid_G1, (count_G1/count_all)*kde_G1, 'r-'); % plot kde for G1 gene angles in red
plot(xgrid_S, (count_S/count_all)*kde_S, 'g-'); % plot kde for S gene angles in green
plot(xgrid_G2, (count_G2/count_all)*kde_G2, 'b-'); % plot kde for S/G2 gene angles in blue
plot(xgrid_M, (count_M/count_all)*kde_M, 'y-'); % plot kde for G2M gene angles in yellow
plot(xgrid_MG1, (count_MG1/count_all)*kde_MG1, 'c-'); % plot kde for MG1 gene angles in cyan
hold off;
legend('top 200 gene angles', 'G1', 'S', 'G2', 'M','M/G1');
axis([0, 2*pi, 0, 0.4]);
% title('kernel density estimation for all top gene angles and each of five subpopulations (us)');
% xlabel('gene angle');
% ylabel('kde');
set(gca,'XTickLabelMode','manual') ;
set(gca,'YTickLabelMode','manual') ;


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