function [prop_P,prop_name]=propAnovaPvalues(Sgrouped, level, levelVar, threshPlotP);
% propertiesAnovaPvalues is used to find significant differences between N groups 
% by univariate one-way unbalanced ANOVA, and show the results by box plots.
% A non-parametric test, Kruskal-Wallis, may also be used
% Program by Steinar Thorvaldsen, steinart@math.uit.no, Nov. 2004. 
% Ref.: http://www.math.uit.no/bi/deltaprot/
% Last changes 27. April 2005.

% Use:     [prop_P,prop_name]=propAnovaPvalues(Sgrouped, level, levelVar, 0.05);

% Input
% Sgrouped: Cell array containing the gene sequence from N
%           different groups of organisms.
%           All seq must be of same length.
% level:    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.
% threshPlotP: Threshold p-value for performing plots. 
% A number (Nprop) of chemical property are read from a property table for aa.

% Output
% prop_P:    Nprop x layers array. Statistical p-values for the observed seq difference.
%            One extra layer is for the total sequence
% prop_name: String array with name of properties

% Labeling 5 groups in 15-character array (****may be changed****):
labels=['Mesophile      '; 'Psychrototerant'; 'Psychrophile   '; 'Group 4        '; 'Group 5        '];

% Reconstructing input data into standard arrays and vectors:
Sint=cell2mat(Sgrouped);
[I,J] = size(Sint); %J must be same for all sequence groups
Ngroups = size(Sgrouped,1); % the number of groups
group=[];
for n=1:Ngroups
    for m=1:size(Sgrouped{n},1); % the number of sequences in group n
        group=[group; labels(n,:)]; %assigning groups lable name in vector
    end
end

% The function is made flexible conserning input parameters:
if nargin < 4
    threshPlotP=0; % default no plotting
end
if nargin < 3
    % Conserved sites will be included in statistical analysis
    levelVar=2*ones(J,1,'int8'); % no variation levels
end
if nargin < 2
    level=ones(J,1,'int8'); % no structural levels
end

levelVar1=find(levelVar<=1); % the conserved sites

D3=max(level); % the number of layers
if D3 > 1 % Constructe D3+1 cells with index vectors for each level:
    D3=D3+1; %append extra space for a total sequence layer
    for n=1:D3-1
        level_index{n} = {find(level==n)};
    end
    level_index{D3} = {find(level>-1)}; %should contain all positions in seq
end 

% Adding up the non-conserved (independent) part of the alignment into 
% a Ix20xD3 matrix for the statistical analysis:
Saa=zeros(I,20,D3); %will contain the number of each aa in each seq and layer
for seqNo=1:I 
    for aaNo=1:20
        current_aa_ind=find(Sint(seqNo,:)==aaNo);
        current_ind=setdiff(current_aa_ind,levelVar1); % fully conserved sites should be excluded
        for layer=1:D3
            layer_ind=cell2mat(level_index{layer});      
            current_ind_layer=intersect(current_ind,layer_ind);
            Saa(seqNo,aaNo,layer)=length(current_ind_layer);
        end
    end       
end

% Read tabular values for 100+ amino acid properties from Excel file:
[prop_v,prop_name]=xlsread('aa_prop60.xls');
[Ip,Jp]=size(prop_v);
prop_name=prop_name(3:Ip+2,22); % to extract names from column nr 22 (top margin=2)
prop_P = zeros (Ip,D3); % the p-values for Ip properties separated in D3 layers

% Perform statistical testing:
hit=[3 5];%[75 79 88 89 90 91 95 109]; %make a list of importante hits after prescanning
for ip=1:Ip %=hit scanning through aa properties
    prop_v_use=prop_v(ip,1:20); %property 21(gaps) not used here
    for layer=1:D3 % Each layer is analysed separately:
        y_prop=zeros(I,1);
        for seqNo=1:I
            % Use element-by-element product, and then compute the mean
            y_prop(seqNo)=sum(prop_v_use.*Saa(seqNo,:,layer))/sum(Saa(seqNo,:,layer)); 
        end
        prop_P(ip,layer)= anova1(y_prop,group,'off'); % perform ANOVA test
        %prop_P(ip,layer)= kruskalwallis(y_prop,group,'off'); % non-parametric test
        if layer ~= D3-1999 & prop_P(ip,layer) < threshPlotP %-999 %0.01 % candidate for an interesing result
            figure
            prop_P(ip,layer)= anova1(y_prop,group);
            title(['Variation in prop ', num2str(ip),' between groups, level ', num2str(layer),', p=',num2str(prop_P(ip,layer))]);
            hold off
        end
    end
end
