function [b,prop_P,prop_name]=propRegressKendall(Sgrouped, level, siteVar, meanBlocks)
% propertiesRegressKendall is used to detect trends in a set of SEVERAL proteins (Nprot) 
% by non-parametric linear regression.
% The observations must be independent, the site must be identically distributed, 
% and mean values for each protein must have symmetric distributions about their means.
% Program by Steinar Thorvaldsen, steinar.thorvaldsen@uit.no, Aug 2008. 
% Ref.: The DeltaProt toolbox at http://www.math.uit.no/bi/deltaprot/
% Last changes 30. Aug 2008.

% Use:     [b,prop_P,prop_name]=propRegressKendall(Sgrouped, level, siteVar, meanBlocks);

% Input
% Sgrouped: Cell array containing the aligned gene sequence from 2 or more
%           different groups of organisms.
%           All sequences must be of same length.
% level:    Vector containin information on 2D or 3D about the aligned sequences. 
%           The number of layers can be 1.
% siteVar: Variation level at each site in the sequence. Conserved sites has siteVar=1.
% meanBlocks: Optional input that may be used to apply test on MEAN values (prefered)
% in blocked strata (=each protein).
% A number,Nprop, of chemical properties to scan. Are read from a property table for aa.

% Output
% b: Nprop x layers x Nprot x 2 array. Regression coefficients in the model Y = b1 + b2 * X + e
% prop_P:    Nprop x layers array. Statistical p-values for the observed seq difference.
% prop_name: String array with name of properties

% Reconstructing the dimensions of the data:
Sint=cell2mat(Sgrouped);
[I,J] = size(Sint); %J must be same for all sequence groups
Spop=ones(I,1); %Vector to store population membership of each sequence
Ngroups = size(Sgrouped,1); % the number of populations
i=0;
for n=1:Ngroups
    for m=1:size(Sgrouped{n},1); % the number of sequences in group n
        i=i+1;
        Spop(i)=n; %assigning population membership in the vector
    end
end

% The function is made flexible conserning input parameters:
if nargin < 2 | isempty(level)
    level=ones(J,1,'int8'); % no structural layers
end
if nargin < 3 | isempty(siteVar)
    siteVar=2*ones(J,1,'int8'); % no site variation levels
end

% Some seqence sites may be excluded from the statistical tests:
siteVar1=find(siteVar<=1); % the conserved sites
%gaps=find(max(cell2mat(Sgrouped))>=21); %the positions with gaps (no 21)
%excl=union(siteVar1,gaps);

D3=max(level); % the number of layers in the level filter
if D3 > 1 % Constructe D3+1 cells with index vectors for each level:
    D3=D3+1; %append extra space for the sum of all other layers
    for n=1:D3-1
        level_index{n} = {find(level==n)};
    end
    level_index{D3} = {(1:J)}; %containing all levels (should gaps be excluded?)
    %level_index{D3} = {find(level<=D3-2)}; %containing all levels of template, gaps excluded   
end 

% Read tabular quantitative values for 60+ amino acid properties from Excel file:
[prop_q,prop_name]=xlsread('aa_prop60.xls');
[Mp,Np]=size(prop_q);
if nargin == 4 & ~isempty(meanBlocks)
    % Performing analysis on the mean value of each protein sequence:
    prop_q(:,21)=NaN; % all gaps (no 21) will then later be excluded
    Nprot=length(meanBlocks); %the number of proteins (blocks)
    Jblocks=sum(meanBlocks);
    if J ~= Jblocks
         error('Length of input arguments Sgrouped and meanBlocks do not match')
    end
    b = zeros(Mp,D3,Nprot,2); % Two coefficients in each regression model (intercept and slope)
else
    b = zeros(Mp,D3,1,2); % less space for coefficients needed in the pooled regression model!
end
prop_name=prop_name(3:Mp+2,22); % to extract names from column nr 22 (top margin=2)
Sprop_q = zeros (I,J);
prop_P = ones(Mp,D3)*NaN; % the P-values for Mp properties separated in D3 layers

important_hit = [4 5 7 8 30 31 36 47 54 55 56 57 58 59 85 86]; %make a list after prescanning
for mp=important_hit %1:Mp %scanning through all properties
    disp (['Testing property: ', num2str(mp),'... '])
    % Substituting property values in all sequences:
    for i=1:I
      Sprop_q(i,:) = prop_q(mp,Sint(i,:)); %property-values stored in each seq position
      % NB: only aa-numbers between 1..21 must existe in sequence!
    end

    % Performing analysis on the MEAN VALUE of each sequence, not on all sites in alignment:
    % (this is the prefered approach)
    if nargin == 4 & ~isempty(meanBlocks)
      % Each layer is analysed separately:
      for layer=1:D3 %=[1 2 3 5]
        Sprop_q_current=ones(I,J)*NaN;
        Sprop_q_mean=ones(I,Nprot)*NaN;
        current_ind=cell2mat(level_index{layer});
        %To exclude fully conserved sites from test:
        %current_ind=setdiff(current_ind,siteVar1);
        Sprop_q_current(:,current_ind)=Sprop_q(:,current_ind);
        start=1;
        stop=0;
        for protein =1:Nprot
          stop=stop+meanBlocks(protein);
          Sprop_q_mean(:,protein)=nanmean(Sprop_q_current(:,start:stop),2);
          X = [ones(size(Spop)) Spop]; % column of 1 because model contains a constant term
          %The backslash operator enables ordinary least sqares estimation:
          b(mp,layer,protein,:) = X\Sprop_q_mean(:,protein); 
          start=stop+1;
        end
        slopes = squeeze(b(mp,layer,:,2)); %converts matrix2vector
        
        % Her we implement 2 alternative methods: Wilcoxon and Kendall (prefered):       
        % The Wilcoxon signed rank test for slope H: b2 <> 0:
        %nodata=ones(length(slopes),1)*NaN;
        %if not(isequalwithequalnans(slopes,nodata)) %to avoid error in signrank
        %  prop_P(mp,layer)= signrank(slopes);
        %end
        
        % As an alternative the Kendall rank test may be used:
        prop_P(mp,layer)= cumKendallTest(Sprop_q_mean,Spop);
        if nanmean(slopes)<0
           prop_P(mp,layer)=-prop_P(mp,layer); %just to indicate direction of slope
        end
      end
    
% -----Here is our old approach where gaps are problematic--------------------------
    else %Performing analysis based on EACH SITE of the alignment.
      %First we compute all the possible number of slopes at a sequence position:
      Suni=union(Spop,Spop); %The unique populations
      SlopesN=0;
      for ii=1:length(Suni)
          biN=length(find(Spop, Suni(ii)));
          SlopesN=SlopesN+biN*(I-biN);
      end
      SlopesN=SlopesN/2; %each slope has been counted twice 
      slopesAll=ones(SlopesN,J)*NaN; %all possible pairs of slopes at pos J from data
      %Next we compute all the possible slopes at each sequence position:
      for j=1:J
        k=0;
        for ii=1:I-1
          for ii2=ii+1:I
              if Spop(ii)~= Spop(ii2)
                k=k+1;
                change_q = Sprop_q(ii2,j)-Sprop_q(ii,j);
                if change_q ~= 0 % to avoid conserved amino acids
                   slopesAll(k,j)=(Sprop_q(ii2,j)-Sprop_q(ii,j))/(Spop(ii2)-Spop(ii));
                end
              end
          end
        end
      end
      % Then each layer is analysed separately:
      for layer=1:D3 %=[1 2 3 5] 
        current_ind=cell2mat(level_index{layer});
        %To exclude fully conserved sites:
        current_ind=setdiff(current_ind,siteVar1);
        siz=prod(size(slopesAll(:,current_ind)));
        %The pooled Theil-Sen slope estimator: 
        b(mp,layer,1,2) = nanmedian(reshape(slopesAll(:,current_ind),siz,1)); % often =0!
        %b(mp,layer,1,2) = nanmean(reshape(slopesAll(:,current_ind),siz,1)); %must be vector
        b(mp,layer,1,1) = NaN; %not relevant output in this case
        % The cumulative Kendall trend test:
        prop_P(mp,layer)= cumKendallTest(Sprop_q(:,current_ind),Spop);        
      end
    end
end
