function [b,prop_P,prop_name]=propRegress(S, X, level, siteVar)
% propertiesRegress is used to detect trends in ONE protein 
% by parametric linear regression. 
% Several taxionomic strata may be included in the model using dummy varables for strata
% Since the sequences are of approximately the same length, the mean propety value 
% for each sequence will have the same normal distribution.
% Program by Steinar Thorvaldsen, steinar.thorvaldsen@uit.no, Aug 2008. 
% Ref.: http://www.math.uit.no/bi/deltaprot/
% Last changes 12. Sept 2008.

% Use:     [b,prop_P,prop_name]=propRegress(S, X, level, siteVar);

% Input
% S:        M x N matrix of observations (the alignment with integers).
% X:        M x r matrix of regressors
% level:    Vector containin information on 2D or 3D about the aligned sequences. 
%           The number of levels can be 1.
% siteVar: Variation level at each site in the sequence. Conserved sites has siteVar=1.

% A number,Nprop, of chemical properties to scan. Are read from a property table for aa.

% Output
% b: Nprop x levels x r matrix of regression coefficients 
% in the model, for example: Y = b0 + b1 * X1 + b2 * I1 +b3 * I2 + e
% prop_P:    Nprop x layers matrix. Statistical p-values for the observed seq difference.
% prop_name: String array with name of properties

% The dimensions of the data:
[M,N] = size(S); %N must be same for all sequences

% The function is made flexible conserning input parameters:
if nargin < 4 || isempty(siteVar)
    siteVar=2*ones(N,1,'int8'); % no site variation levels
end
if nargin < 3 || isempty(level)
    level=ones(N,1,'int8'); % no structural layers
end
[M1,r] = size(X);
if M ~= M1
    error('propRegress: Matrix with observations must match the matrix of regressors');
end
X = [ones(M1,1) X]; % adds column of 1 because model contains a constant term***

% Seqence sites to 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 2D/3D 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:N)}; %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');
prop_q(:,21)=NaN; % all gaps (no 21) will then later be excluded
[Mp,Np]=size(prop_q);
b = zeros(Mp,D3,r+1); % coefficients in each regression model (b0=intercept and b1=slope)
prop_name=prop_name(3:Mp+2,22); % to extract names from column nr 22 (top margin=2)
Sprop_q = ones (M,N)*NaN;
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]; %make a list after prescanning
important_hit =[55 56 54 5 57 4 58 59 30 7 31 8 47 36 85]; %same order as in table 2, 3
%important_hit =[15 16 18 85];
for mp=important_hit %1:Mp %scanning through all properties
    disp (['Testing property: ', num2str(mp),'... '])
    % Substituting property values in all sequences:
    for m=1:M
      Sprop_q(m,:) = prop_q(mp,S(m,:)); %property-values stored in each seq position
      % NB: only aa-numbers between 1..21 must existe in sequence!
    end
    % Each layer is analysed separately:
    for layer=1:D3 %=[1 2 3 5]
        Sprop_q_current=ones(M,N)*NaN;
        Sprop_q_mean=ones(M,1)*NaN;
        current_ind=cell2mat(level_index{layer});
        %To eventually exclude fully conserved sites from test:
        %current_ind=setdiff(current_ind,siteVar1);
        Sprop_q_current(:,current_ind)=Sprop_q(:,current_ind);
        Sprop_q_mean=nanmean(Sprop_q_current,2);
        %Ordinary least squares estimation (Matlab function regstats may also be used):
        [b(mp,layer,:),bint,res,resint,stats] = regress (Sprop_q_mean,X);
        % We must compute P-value for slope b1=0 from 95% confidence interval(bint):
        % M is the sample size, and regression prediction is t-distributed
        % with M-2 degress of freedom:
        SDb1=abs(bint(2,1)-b(mp,layer,2))/tinv(0.975,M-2);
        t=abs(b(mp,layer,2)/SDb1); % use abs for two sided P-value
        prop_P(mp,layer)=2*(1-tcdf(t,M-2));
        if b(mp,layer,2)<0 % the slope
           prop_P(mp,layer)=-prop_P(mp,layer); %just to indicate direction of slope
        end        
        %As an alternative the Kendall rank test may be used:
        %prop_P(mp,layer)= cumKendallTest(Sprop_q_mean,X(:,2));
    end
end