function [] = substPairsPlot(subst, p, alpha);
% substPairsPlot is used to plot the substitution matrix of amino acids,
% gaps (insertions and deletions) may be included.
% Program by Steinar Thorvaldsen, april 2004. Last changes 25. Nov 2005.
% Remark: This script works in Matlab v6.5 and 7.1, but not in v7.01

% Use:   substPairsPlott(subst,p,alpha);

% Input
% subst: 20x20 (or 21x21) matrix containing the aa subst (and indels)
% p    : 20x20 (or 21x21) matrix containing statistical p-values.
%        The sign of p only used to indicate direction of change, 
%        p positive: over representation, p negative: under representation.
% alpha: also performs special tagging at the significance level p, for
% example alpha = 0.01

% Maximum disk size to use in plot, may be changed:
big=300; %Biggest circle area in grid plott.

if nargin < 3
    alpha=0.01; % default value
end

if nargin > 1
    if size(subst)~=size(p)
        error('Substitution matrix and p-value matrix must be of same size.')
    end;
end

[A,B,D3] = size(subst);
if ((A ~= 20)|(B ~= 20)) & ((A ~= 21)|(B ~= 21)) 
   error('Substitution matrix must be of size 20x20 or 21x21(with indels).')
end;
if A == 20 %=B
    aa_letter=char('A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V','');
else % A=B=21(with gaps: ~)
    aa_letter=char('A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V','~','');
end

% Interchanging the axis to make the plot more conventional: FROM x-axis TO y-axis.
for d=1:D3
    subst(:,:,d)=subst(:,:,d)';
    p(:,:,d)=p(:,:,d)';
end

A2 =A*A;
% Each layer is analysed separately:
for d=1:D3
    figure %image, pcolor or scatter functions may all be used here.
    %imagesc(subst(:,:,d)-diag(diag(subst(:,:,d))));%diagonal elements set to zero
    %To plot by scatter, we must first convert the input 20x20 subst matrix to vector form:
       gridpoint=zeros(A);
       for i=1:A gridpoint(i,:)=1:A; end
       x=zeros(A2,1); y=zeros(A2,1); z=zeros(A2,1);
       x(:)=gridpoint(:,:);
       y(:)=gridpoint(:,:)'; %flipud change direction of y-axis for plot
       z(:)=subst(:,:,d)-diag(diag(subst(:,:,d)));%diagonal elements set to zero
    scale=big/max(z); %may be changed to have disks of desired size in plott
    scatter(x(:),y(:),(z(:)+eps)*scale,(z(:)+eps)*scale,'o','filled');
    
    if nargin > 1 % statistical significance are also to be marked at the plott
        hold on;
        % The diagonal elements (non-substitutions) are not to be considered:
        p_matrix=p(:,:,d)-diag(diag(p(:,:,d)))+diag(diag(ones(A)));%all diagonal elements set to 1
        p_lowIndex=find(abs(p_matrix) < alpha); %returns the linear indices corresponding to prob<alpha
        %subst_diff= subst(:,:,d)-subst(:,:,d)';  %the changes in numbers
        upIndex=find(p_matrix>=0);
        downIndex=find(p_matrix<0);
        %The significant favoured substitutions are marked with upward-pointing triangles:
        currentIndex=intersect(upIndex,p_lowIndex);
        scatter(x(currentIndex),y(currentIndex),big,'^k'); 
        %The significant downward substitutions are marked with downward-pointing triangles:
        currentIndex=intersect(downIndex,p_lowIndex);
        scatter(x(currentIndex),y(currentIndex),big,'vk');
    end
    
    title(['Substitutions from mesophiles --> psychrophiles, layer ', num2str(d)]);
    xlabel('Replaced amino acid');
    ylabel('Replacement amino acid');
    axis([0 A+1 0 A+1]);
    box on;
    set(gca,'XTick',1:1:A);
    set(gca,'XTickLabel',aa_letter); %labeling plot in order of the aa
    set(gca,'YTick',1:1:A);
    set(gca,'YTickLabel',aa_letter);
    
    % Becauce of the scaling above, the axis of the colorbar must be adjusted to
    % the right scale:
    h_ax=colorbar;
    old_labels = get(h_ax,'YTicklabel'); % current labels
    [nr,nc] = size(old_labels);
    new_labels = cell(nr,1); %convert from strings to numbers
    old_cell = cellstr(old_labels);
    for j=1:nr %re-scale labels
        new_labels{j}=num2str(round(str2num(old_cell{j})/scale)); %back to numbers
    end 
    set(h_ax,'YTicklabel',new_labels);
    hold off;
end