function [] = substPvalueSort(subst, p, p2, topN);
% substPvalueSort is used to sort the substitution matrix of amino acids,
% in order of increasing P-value in the matrix p.
% Program by Steinar Thorvaldsen, april 2005. Last changes 16. Jan 2006.

% Use:   substPvalueSort(subst,p,p2,topN);

% 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 is only used to indicate direction of change, 
%        p positive: over-representation, p negative: under-representation.
% p2   : optional P-values from different background test
% topN : only repport top N number of counts

if nargin < 2 
   error('Not enought input data.')
elseif nargin < 3 
    if size(subst)~=size(p)
      error('Substitution matrix and p-value matrices must be of same size.')
    end
elseif nargin < 4
    if size(subst)~=size(p2)
      error('Substitution matrix and p-value matrices must be of same size.')
    end;
end;
if nargin < 4
    topN=20; % default value
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

A2=A*A;
if topN > A2
    topN=A2-A; %not including diagonal elements
end
% Each layer is analysed separately:
for d=1:D3
    p_matrix=zeros(A,B);
    % 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
    q = reshape(p_matrix,A2,1); %reshape matrix to a vector
    [q,IX]=sort(abs(q)); %IX is the linear indices corresponding to the sorted matrix
    disp (['The most biased substitution pairs (SP) at level: ',num2str(d)])
    disp ('SP   Forward Reverse P_values (with +/- direction)')
    if nargin > 2 % two sets of statistical P-values are to be reported
        p_matrix2=p2(:,:,d)-diag(diag(p2(:,:,d)))+diag(diag(ones(A)));
    end
    for N=1:topN
     [i j] = ind2sub([A B], IX(N)); %convert from linear indexing to row-column subscripts
     if nargin > 2 % two sets of statistical P-values are to be reported
       fprintf('%c%c: %6.0f %6.0f    P= %9.7f (%9.7f)\n',aa_letter(i),aa_letter(j),subst(i,j,d),subst(j,i,d),p_matrix(i,j),p_matrix2(i,j))
       %disp ([aa_letter(i),aa_letter(j),':   ', num2str(subst(i,j,d)),'     ',num2str(subst(j,i,d)),'       P= ',num2str(p_matrix(i,j)),' (',num2str(p_matrix2(i,j)),')'])    
     else
       fprintf('%c%c: %6.0f %6.0f    P= %9.7f\n',aa_letter(i),aa_letter(j),subst(i,j,d),subst(j,i,d),p_matrix(i,j))
       %disp ([aa_letter(i),aa_letter(j),':   ', num2str(subst(i,j,d)),'     ',num2str(subst(j,i,d)),'       P= ',num2str(p_matrix(i,j))])      
     end
    end
    disp (' ') %empty line
end