function substBin = substPairs2bin(substPairs, method);
% substPairs2bin is used to reduces the amino acid substitution matrix to 
% fewer categories based on selected grouping of related amino acids
% Program by Steinar Thorvaldsen, steinart@math.uit.no, July 2006. 
% Ref.: http://www.math.uit.no/bi/deltaprot/
% Last changes 22. Dec 2006.

% Input:
% substPairs: integer array containing the number of each aa substitution.
% method: defining the new grouping of amino acids:
% 1:hydrophobic, 2:charged, 3:positive, 4:negative, 5:polar,
% 6:aromatic, 7:aliphatic, 8:tiny, 9:small, ...
% Output:
% subst2Bin:  2x2 matrix with the accumulated bin counts.
% The following mapping protocol is used:
% A R N D C Q E G H I  L  K  M  F  P  S  T  W  Y  V  - 
% 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
% where - are gaps.
% Use:
% substCharged = substPairs2bin(substPairs, 2);

% The function is made flexible conserning input parameters:
if nargin < 2
    method=1;
end

[M,N,P] = size (substPairs);
if (M ~= 20 & M ~= 21) | (N ~= 20 & N ~= 21)
    error('substPairs2bin: Substitution matrix must be of size 20x20 or 21x21');
end

method = round(method);
switch method
    case 1 % hydrophobic
        bin1=[1 8 10 11 13 14 15 18 20];
    case 2 % charged
        bin1=[2 4 7 9 12];
    case 3 % positive
        bin1=[2 9 12]; %9(=H) may be removed...
    case 4 % negative
        bin1=[4 7];
    case 5 % polar
        bin1=[3 5 6 9 16 17 18 19];
    case 6 % aromatic
        bin1=[9 14 18 19];
    case 7 % aliphatic
        bin1=[10 11 20];
    case 8 % tiny
        bin1=[1 8 16];
    case 9 % small
        bin1=[1 3 4 5 8 15 16 17 20];
    case 10 % pos or polar
        bin1=[2 3 5 6 9 12 16 17 18 19];
    case 11 % neg or polar
        bin1=[3 4 5 6 7 9 16 17 18 19];
    otherwise % No defined method
        error('substPairs2bin:UnknownMethod', ...
          'This ''method'' parameter is not known.');    
end
bin2=setdiff([1:M],bin1); %the rest of the amino acids

substBin=zeros(2,2,P);
for p=1:P
    %removing diagonal elements (non-substitutions):
    sP=substPairs(:,:,p)-diag(diag(substPairs(:,:,p))); 
    substBin(1,1,p)=sum(sum(sP(bin1,bin1)));
    substBin(2,1,p)=sum(sum(sP(bin2,bin1)));
    substBin(1,2,p)=sum(sum(sP(bin1,bin2)));
    substBin(2,2,p)=sum(sum(sP(bin2,bin2)));
    %For alternative pooling into 1x2 matrices:
    %substBin(1,1,p)=sum(sum(sP(bin1,bin2))); %subst bin1->bin2
    %substBin(1,2,p)=sum(sum(sP(bin2,bin1))); %subst bin2->bin1
    %substBin(1,1,p)=sum(sum(sP(:,bin1))); %subst x->bin1
    %substBin(1,2,p)=sum(sum(sP(:,bin2))); %subst x->not(bin1)
end
