% DeltaProt: protein comparision, main program script. 
% A comparative analysis of aligned, homolog protein sequences.
% Program by Steinar Thorvaldsen, steinar.thorvaldsen@uit.no, April. 2009.
% Ref.: S.Thorvaldsen and E.Ytterstad: "Environmental adaptation of proteins: 
% Regression models with simple physicochemical properties".
% Computational Biology and Chemistry. Vol. 33 (5) 2009, Pages 351-356.
% Last uppdated 8.Oct 2009

% Input:
% Fasta files with m aligned protein sequences extracted from different organisms.
% (fasta format). The files must be stored in a subdirectory like \memb\
% To each organisme biological and phenotype parameters may be attached 
% to the header of the sequence.

% Auxilary predictions to be runned beforehand (external servers)
% SABLE:    Relative solvent accessibility of the residues
%           (http://sable.cchmc.org/)
%           (0-completely buried, 9-fully exposed) 
%           Confidence level of the RSA prediction (0-low confidence, 9-highly confident) 

clear all;
% You may use the global constant 'maxlevel' to force the analyses of all
% genes to be done with the same number of 2D and 3D levels:
% (not relevant when only one alignment is to be analysed)
global maxlevel;
    maxlevel=4; % levels of 2D/3D in analyses
    
%---------------------- Read and convert data --------------------------
aa_seq = dir('./memb/*.stn'); % reads the filname from data directory into a structure array with fields
maxNo=size(aa_seq,1); % the number of files in this particular case
M=12; %the number of data lines represented in this data set (the fasta-files)

% Defining sequence populations of all 7 sequences (grouped by temperature):
maxPop=3;
t31_40=[2 3 4]; % the meophile sequences, nr 1 and 2 are a copies
t21_30=[5]; % the intermediate sequence
t1_20=[6 7]; % the psycrophile sequences.
extra_lines=5; %line 8..12 reserved for 2D/3D information and optional template seqence

%Here dim 2 and 4 must match with data! 
aaDeltafreq2D = zeros (maxNo,maxPop-1,20,maxlevel+1); % for storing changes in frequencies
aaDeltafreq3D = zeros (maxNo,maxPop-1,20,maxlevel+1); % for storing changes in frequencies
proteinBlock = zeros(maxNo,1); % for storing the length of each protein sequence block
Stotal=[];
levelVartotal=[];
level2Dtotal=[];
level3Dtotal=[];

for fileNo=1:maxNo
    % L is the number of lines in the alignment file, n is the length, 
    % S is the character sequence array, and names are the headers:
    [L, n, S, names] = readfasta(strcat('./memb/',aa_seq(fileNo).name)); %strcat adds the two strings
    proteinBlock(fileNo) = n; % length of this alignment
    Sint=aa2int(S(1:L-extra_lines,:)); % integer number representation of alignment S
    levelVar = findSeqSiteVar(Sint); % finding the variation and conserved sites (Var=1) along ALL sequences
    % Computational filters added at the top or bottom of the alignment:
    level2D=dssp2int(S(9,:)); % secondary structure(1,2,3)
    level3D=wasa2int(S(11,:)); % internal/external amino acids(only levels 1,2,3)
    Stotal=[Stotal,Sint]; % store data by just extending the arrays in dim 2
    levelVartotal=[levelVartotal,levelVar];
    level2Dtotal=[level2Dtotal,level2D]; 
    level3Dtotal=[level3Dtotal,level3D];  
end

% We use cell arrays to store dissimilar size of data:
Sgrouped= {Stotal(t31_40,:); Stotal(t21_30,:); Stotal(t1_20,:)};

% 3D structure of protein (exsample)
substMP3D = substPairsCount(Stotal(t31_40,:), Stotal(t1_20,:), level3Dtotal, 1, 1);
substMM3D = substPairsCount(Stotal(t31_40,:), Stotal(t31_40,:), level3Dtotal, 1, 0);

PvalueMPvsMM3D = substPairsPvalues(substMP3D,substMM3D);
substPairsPlot(substMP3D,PvalueMPvsMM3D, 0.01);
substPvalueSort(substMP3D,PvalueMPvsPM3D);
substCountSort(substMP3D,PvalueMPvsPM3D);
FDR (PvalueMPvsMM3D, 0.05)

substMP3D = substPairsCount(Stotal(t31_40,:), Stotal(t1_20,:), level3Dtotal, 1, 0);
substMM3D = substPairsCount(Stotal(t31_40,:), Stotal(t31_40,:), level3Dtotal, 1, 0);
PvalueMPvsMM3D = substPairsPvalues(substMP3D,substMM3D);
substPairsPlot(substMP3D,PvalueMPvsMM3D, 0.01);
FDR (PvalueMPvsMM3D(:,:,5), 0.05)
%-----------Regression analyses of property trends --------------
[b2D,regressP2D,prop_name]=propRegressKendall(Sgrouped, level2Dtotal, levelVartotal, proteinBlock);
[b3D,regressP3D,prop_name]=propRegressKendall(Sgrouped, level3Dtotal, levelVartotal, proteinBlock);
%[b3,regressSiteP3D,prop_name]=propRegressKendall(Sgrouped, level3Dtotal, levelVartotal);
yer = [4 5 7 8 30 31 36 47 54 55 56 57 58 59 85];
yertab = [55 56 54 5 57 4 58 59 30 7 31 8 47 36 85];
xer = [1 2 3 5];
FDR (abs(regressP3D(yertab,xer)), 0.05) %false discovery rate

% For ploting slopes of one property:
Property47level5_Allprotein_Slope=sort(squeeze(b3D(47,5,:,2))); %for ploting
figure;
bar (Property47level5_Allprotein_Slope); % Fig. 1 in Reference article