function [aaComp, aaStdGroup1, aa_total] = aaFreq(Sgrouped, level, levelVar);
%function [aaComp, aaStdGroup1, aa_total, aa_freq] = aaFreq(Sgrouped, level, levelVar);
% aaFreq is used to find the frequences of amino acids, conserved sites may be excluded
% Program by Steinar Thorvaldsen, april 2004. Last changes 22. Dec 2006.

% Use:      
%           [aaComposition, aaStd1] = aaFreq(Sgrouped, level, levelVar)
% Input
% Sgrouped: cell array containing the gene sequence from 2 or more different groups of organisms
% level:    a vector containin information on 2D or 3D about the aligned sequences. 
%           The number of layers may be 1 as for unaligned sequences.
% levelVar: Amino acid variations at each site in the alignment.
%           The variation =1 at conserved sites, and they may be excluded to
%           avoide statistical dependence.
% Output
% aaComp:   Ngroups x 20 x levels, a level-separated matrix with aa countings from the N groups.
%           Gaps are neglected.
% aaStdGroup1: Standard deviation of group 1 (the baseline group).
% aa_total: M x 20 x levels with counts of aa, where M is the number of seq.
%           The seq in group 1 will come first, then group 2 etc.
% REMARK:   all output data are sorted by group!

% Reconstructing the dimension of the data:
Ngroups = size(Sgrouped,1); % the number of groupes
for n=1:Ngroups
    nGroup(n)=size(Sgrouped{n},1); % the number of sequences in each group
end
M = sum(nGroup); % the total number of sequences
[M1,N1] = size(cell2mat(Sgrouped(1))); %N1 must be same in all seq groups

% The function is made flexible conserning input parameters:
if nargin < 3
    % Conserved sites will be included in statistical analysis
    levelVar=2*ones(N1,1,'int8'); % no variation levels
end
if nargin < 2
    level=ones(N1,1,'int8'); % no structural levels
end
D3=max(level); % the number of layers in the level filter

levelVar1=find(levelVar<=1); % the conserved sites
site_ind=[1:N1];
current_ind=setdiff(site_ind,levelVar1); % fully conserved sites should be excluded
S=cell2mat(Sgrouped);

% Starting to compute frequences of amino acids:
aa_total = aaCount(S(:,current_ind), level(current_ind)); % aa compositions for each seq from function aaCount
if D3 > 1
    aa_total=cat(3,aa_total,sum(aa_total,3)); %Concatenate two arrays by appending sums as extra layer
    D3=D3+1;
end
aaComp = zeros (Ngroups,20,D3); % the frequency data for N groups separaten in D3 layers
aaStdGroup1 = zeros (20,D3); % standard deviation for group 1

% Each layer is analysed separately:
aa_freq = zeros (M,20,D3);
for layer=1:D3
    % Transformeing the matrix to percentage:
    for m=1:M
        aa_sum = sum(aa_total(m,:,layer),2); %the total number of aa
        if aa_sum > 0
            % here we use elementwise division:
            aa_freq(m,:,layer) = 100*aa_total(m,:,layer)./aa_sum;
        end
    end
    aa_freq_group = zeros (Ngroups,20);
    % To finde the mean frequence and standard deviation of the baseline group (1):
    aa_freq_group(1,:)=mean(aa_freq(1:nGroup(1),:,layer),1); %mean values in dim 1
    aa_std_group1=std(aa_freq(1:nGroup(1),:,layer),0,1); %flag=0, dim=1
    % To find the variation in mean aa composition in the other groups of sequences:
    startSeq=nGroup(1);
    for n=2:Ngroups
        stopSeq=startSeq+nGroup(n);
        aa_freq_group(n,:)=mean(aa_freq(startSeq+1:stopSeq,:,layer),1);
        startSeq=stopSeq;
    end
    aaComp(:,:,layer) = aa_freq_group(:,:); % returne data array
    aaStdGroup1(:,layer) = aa_std_group1(:);
    
    % To run a statistical test if the frequences are from the same multinomial distribution:
    % P = chi2Test(aa_total(:,:,layer))
end %layer loop
