
%
% Load the data
%

% expressions of the genes, a 1475x20 matrix
exp = load('NODE_DATA.dat');

% links, a matrix of 7123x2, where each row tells one interaction
links = load('LINKS.dat');

% convert links to a binary matrix of 1475x1475
linkmat = zeros(size(exp,1));
for i=1:size(links,1);
  linkmat(links(i,1),links(i,2)) = 1;
end;

classes = load('GOclasses.txt');


%
% Cluster the expressions using K-means
%

Klist = [5,10,20,40];

% Some matrices to store results
fres = zeros(length(Klist),5);
fres2 = zeros(length(Klist),5);

% Go through different values of K
for k=1:length(Klist);

  K = Klist(k);
  [k,K]
  [idx1,c] = kmeans(exp, K,'distance','sqEuclidean');

  % Visualize the result
  [foo,order1] = sort(idx1);
  imagesc(exp(order1,:)); colormap('gray')
  display('press space to continue')
  pause

  % Look at the first two mean vectors
  plot(c(1,:)); hold on; plot(c(2,:),'r');
  pause
  

  [idx2,c] = kmeans(exp, K,'distance','correlation');

  %
  % Compare to GO classes
  %
  res1 = compareToGO(idx1,classes);
  res2 = compareToGO(idx2,classes);



  %
  % Hierarchical clustering
  %

  dmat = pdist(exp,'correlation');
  Z = linkage(dmat,'complete');
  T = cluster(Z,'MaxClust',K);

  % does not work for large clusterings
  %dendrogram(Z,30);

  %
  % Add some link information
  %

  % Combination by multiplying disctances with constants
  % below 1 if there is short path through links
  smoothed = smoothLinks(linkmat,3);
  combined = ones(size(smoothed));
  combined(find(smoothed==1)) = 0.05;
  combined(find(smoothed==2)) = 0.15;
  combined(find(smoothed==3)) = 0.35;
  for i=1:size(combined,1);
    combined(i,i)=0;
  end;
  
  dmat2 = squareform(combined,'tovector');
  dmat_both = dmat.*dmat2;
  Zb = linkage(dmat_both,'complete');
  Tb = cluster(Zb,'MaxClust',K);

  %
  % Comparisons
  %
  res3 = compareToGO(T,classes);
  res4 = compareToGO(Tb,classes);

  % Random clustering to act as a baseline
  cl = round(rand(length(T),1)*K+0.5);

  res5 = compareToGO(cl,classes);

  th = 0.0001;

  fres(k,:)=[length(find(res1<th)),length(find(res2<th)),length(find(res3<th)),length(find(res4<th)),length(find(res5<th))]/prod(size(res1));

  fres2(k,:)=[length(find(min(res1')<th)),length(find(min(res2')<th)),length(find(min(res3')<th)),length(find(min(res4')<th)),length(find(min(res5')<th))]/K;

end;

display('Portion of cluster/GOclass pairs that were enriched. First two K-means clusterings with different distance measures then two hierarchical clustering results with and without link information, and finally random baseline.');
display('The rows are solutions with different numbers of clusters, and the columns indicate GO classes.');
fres
display('The same, but telling the portion of clusters that have at least one enriched GO class.');
fres2
