function [] = substCountSort(subst, p, topN);
% substCountSort is used to sort the substitution matrix of amino acids,
% in order of the most frequent counts.
% Program by Steinar Thorvaldsen, april 2005. Last changes 15. Jan 2006.

% Use:   substCountSort(subst,p,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.
% topN:  only repporte top N number of counts

if nargin < 3
    topN=10; % 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

A2=A*A;
% Each layer is analysed separately:
for d=1:D3
    z=zeros(A,B);
    z(:,:)=subst(:,:,d)-diag(diag(subst(:,:,d)));%diagonal elements set to zero
    z = reshape(z,A2,1);
    [z,IX]=sort(abs(z),'descend'); %IX is the linear indices corresponding to the sorted matrix
    disp (['Most frequent substitution pairs (SP) in numbers at level: ',num2str(d)])
    disp ('SP  Forward Reverse  P_value (with +/- direction) ')
    if nargin > 1 % statistical P-values are also to be rapported
        % The diagonal elements (non-substitutions) are not to be considered:
        p_matrix=p(:,:,d)+diag(diag(ones(A)));%all diagonal elements set to 1
        subst_diff= subst(:,:,d)-subst(:,:,d)';  %the changes in numbers
        for N=1:topN
            [i j] = ind2sub([A B], IX(N)); %convert from linear indexing to row-column subscripts
            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
        disp (' ') %empty line
    end
end