%a matlab script
%correlation analysis of AxB RI lines
%cardiovascular traits
%GA Churchill January 2003

%load data
[data, trait, line] = tblread('means2.dat','\t');
ntraits = size(data,2);

%look at correlations using color-scaled matrix plot
cc = corrcoef(data);
figure
matplot(abs(cc))
tick = [1:ntraits]+0.5;
set(gca,'XTick',tick,'YTick',tick,'XTickLabel',trait,'YTickLabel',trait)

%look at pairwise scatterplots
k=1;
for i = 1:ntraits-1
   for j = i+1:ntraits
   	if mod(k,6) == 1
     	 figure
   	end
   	subplot(3,2,mod(k-1,6)+1)
      plot(data(:,i),data(:,j),'.')
      xlabel(trait(i,:))
      ylabel(trait(j,:))
      axis square
      k = k+1;
   end
end

%permutation test for maximum corr coef
[r,c] = size(data);
permdata = data;
nperms = 10000;
maxcc = zeros(nperms,1);
for k = 1:nperms;
	for i = 2:c
   	permdata(:,i) = data(randperm(r),i);
	end
	ccperm = corrcoef(permdata);
   maxcc(k) = max(max(abs(ccperm-eye(ntraits))));
end
maxcc = sort(maxcc);
maxcc(floor([0.900 0.950 0.990]*nperms)')

%clean up work space
close all
clear ans tick r c k tick ccperm permdata i j k 

%cluster the phenotypes 
% use 1-abs(cc) distance and Ward methods to cluster
dist = zeros(1,(ntraits-1)*ntraits/2)
k = 1;
for i = 1:ntraits-1
   for j = i+1:ntraits
      dist(k) = 1-abs(cc(i,j));
      k = k+1;
   end
end
%z = linkage(dist,'ward');
z = linkage(dist);
dendrogram(z)
%tord = [4 5 11 3 10 6 8 13 12 1 2 9 7];
tord = [4 5 3 10 6 11 13 1 2 9 8 7 12];
text([1:ntraits],zeros(1,ntraits),trait(tord,:), ...
    'HorizontalAlignment','right', ...
    'VerticalAlignment','bottom',...
    'Rotatio',-90);

%redraw correlation matrix traits ordered as in dendrogram 
recc = corrcoef(data(:,tord));
figure
matplot(abs(recc))
tick = [1:12]+0.5;
set(gca,'XTick',tick,'YTick',tick,'XTickLabel',trait(tord,:),...
   'YTickLabel',trait(tord,:))

%draw matrix of signficant connections
signif = (abs(recc) > maxcc(0.9*nperms));
signif = signif + (abs(recc) > maxcc(0.95*nperms));
signif = signif + (abs(recc) > maxcc(0.99*nperms));
figure
matplot(signif)
colormap(flipud(gray))
tick = [1:12]+0.5;
set(gca,'XTick',tick,'YTick',tick,'XTickLabel',trait(tord,:),...
   'YTickLabel',trait(tord,:))

clear ans i j k recc tick tord dist z signif
whos

%the traits are
% PWTh  
% SWTh  
% EDD   
% ESD   
% FS        = (EDD-ESD)/EDD  * 100
% LVMASS    = 1.06 * ((EDD+PWTh+SWTh)^3-(EDD)^3)
% BW    
% LV/BW 
% TH/R      =(PWTh + SWTh)/EDD 
% SV        EDD^3-ESD^3  
% HR    
% EST   
% CO        = SV * HR /1000

%permutation test for maximum corr coef using calculated traits
[r,c] = size(data);
permdata = data;
nperms = 100000;
ccperm = zeros(13,13,nperms);
for k = 1:nperms;
    
    %permute primary variables
	for i = [2,3,4,7,8,11,12];
   	permdata(:,i) = data(randperm(r),i);
	end
    
    %calculated vars
    permdata(:,5) = 100*(permdata(:,3)-permdata(:,4))./permdata(:,3); 
    permdata(:,6) = 1.06*((permdata(:,1)+permdata(:,2)+permdata(:,3)).^3-permdata(:,3).^3);  
    permdata(:,9) = (permdata(:,1)+permdata(:,2))./permdata(:,3);  
    permdata(:,10) = permdata(:,3).^3-permdata(:,4).^3;  
    permdata(:,13) = (permdata(:,10).*permdata(:,11))/1000;  

    ccperm(:,:,k) = corrcoef(permdata);
end
ccsort = sort(abs(ccperm),3);

calcdata = data;
calcdata(:,5) = 100*(calcdata(:,3)-calcdata(:,4))./calcdata(:,3); 
calcdata(:,6) = 1.06*((calcdata(:,1)+calcdata(:,2)+calcdata(:,3)).^3-calcdata(:,3).^3);  
calcdata(:,9) = (calcdata(:,1)+calcdata(:,2))./calcdata(:,3);  
calcdata(:,10) = calcdata(:,3).^3-calcdata(:,4).^3;  
calcdata(:,13) = (calcdata(:,10).*calcdata(:,11))/1000; 
cc = corrcoef(calcdata);

% generate significant pairs list at 0.05 bonferroni adjusted
thresh = squeeze(ccsort(:,:,99936));
tmp = ccsort([1,2,3,4,7,8,11,12],[1,2,3,4,7,8,11,12],99936);
mean(tmp(:))
signif = abs(cc) > thresh;
spy(signif)
[i j] = find(signif);
idx = i(find(i<j));
jdx = j(find(i<j));
[trait(idx,:), trait(jdx,:)]


%generate significant pairs list at 0.10 bonferroni
thresh = squeeze(ccsort(:,:,99872));
tmp = ccsort([1,2,3,4,7,8,11,12],[1,2,3,4,7,8,11,12],99872);
mean(tmp(:))
signif = abs(cc) > thresh;
figure
spy(signif)
[i j] = find(signif);
idx = i(find(i<j));
jdx = j(find(i<j));
[trait(idx,:), trait(jdx,:)]

%generate significant pairs list at 0.01 nominal
thresh = squeeze(ccsort(:,:,99000));
tmp = ccsort([1,2,3,4,7,8,11,12],[1,2,3,4,7,8,11,12],99000);
mean(tmp(:))
signif = abs(cc) > thresh;
figure
spy(signif)
[i j] = find(signif);
idx = i(find(i<j));
jdx = j(find(i<j));
[trait(idx,:), trait(jdx,:)]

%generate significant pairs list at 0.05 nominal
thresh = squeeze(ccsort(:,:,95000));
tmp = ccsort([1,2,3,4,7,8,11,12],[1,2,3,4,7,8,11,12],95000);
mean(tmp(:))
signif = abs(cc) > thresh;
figure
spy(signif)
[i j] = find(signif);
idx = i(find(i<j));
jdx = j(find(i<j));
[trait(idx,:), trait(jdx,:)]


%clean up
close all
clear ans tick r c k tick ccperm permdata i j k nperms