
% This code was developed using R2016a Matlab version (windows platform)
% This Program ranks the features k-spaced amino acid pairs which are CKSAAP features of peptides on Xu's dataset using the F-score method
% This was carried out to make sure our code is correct by comparing the rank achieved in our code with what was achieved in CKSAAP_PhoglySite work
% The feature ranks can be obtained from the variable named 'Feature_Rank' after the code has executed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First Phase: Calculate the number of k-spaced amino acid paris in each peptide
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load Xu_Data

Amino_acids = {'A' 'C' 'D' 'E' 'F' 'G' 'H' 'I' 'K' 'L' 'M' 'N' 'P' 'Q' 'R' 'S' 'T' 'V' 'W' 'X' 'Y'};

% k = 0
for h = 1:size(Xu_Data, 1)
    sequence = Xu_Data{h,2};
    feature_num = 0;
    for i = 1:size(Amino_acids, 2)
        for j = 1:size(Amino_acids, 2)
            feature_num = feature_num + 1;
            s = strcat(Amino_acids{i}, Amino_acids{j}); % concatenate 2 amino acids 
            match = strfind(sequence, s); % find the concatenated amino acids in the sequence/peptide
            Num = size(match, 2); % store the number of matches
            Xu_Data{h,feature_num + 4} = Num; % save the number as its feature
            feature_label {1,feature_num + 4} = s;
        end
    end
end


% k = 1
for h = 1:size(Xu_Data, 1)
    sequence = Xu_Data{h,2};
    feature_num = 0;
    for i = 1:size(Amino_acids, 2)
        for j = 1:size(Amino_acids, 2)
            feature_num = feature_num + 1;
            s = strcat(Amino_acids{i}, '(?=.', Amino_acids{j},')'); % concatenate 2 amino acids. Dot is added to indicate there can be any amino acid in between (concatenate to come to this form for e.g. X(?=.X))
            match = regexp(sequence, s); % find the concatenated amino acids in the sequence/peptide
            Num = size(match, 2); % store the number of matches
            Xu_Data{h,feature_num + 445} = Num; % save the number as its feature
            feature_label {1,feature_num + 445} = s;
        end
    end
end

% k = 2
for h = 1:size(Xu_Data, 1)
    sequence = Xu_Data{h,2};
    feature_num = 0;
    for i = 1:size(Amino_acids, 2)
        for j = 1:size(Amino_acids, 2)
            feature_num = feature_num + 1;
            s = strcat(Amino_acids{i}, '(?=..', Amino_acids{j},')'); % concatenate 2 amino acids. Dot is added to indicate there can be any amino acid in between (concatenate to come to this form for e.g. X(?=..X))
            match = regexp(sequence, s); % find the concatenated amino acids in the sequence/peptide
            Num = size(match, 2); % store the number of matches
            Xu_Data{h,feature_num + 886} = Num; % save the number as its feature
            feature_label {1,feature_num + 886} = s;
        end
    end
end

% k = 3
for h = 1:size(Xu_Data, 1)
    sequence = Xu_Data{h,2};
    feature_num = 0;
    for i = 1:size(Amino_acids, 2)
        for j = 1:size(Amino_acids, 2)
            feature_num = feature_num + 1;
            s = strcat(Amino_acids{i}, '(?=...', Amino_acids{j},')'); % concatenate 2 amino acids. Dot is added to indicate there can be any amino acid in between (concatenate to come to this form for e.g. X(?=...X))
            match = regexp(sequence, s); % find the concatenated amino acids in the sequence/peptide
            Num = size(match, 2); % store the number of matches
            Xu_Data{h,feature_num + 1327} = Num; % save the number as its feature
            feature_label {1,feature_num + 1327} = s;
        end
    end
end

% k = 4
for h = 1:size(Xu_Data, 1)
    sequence = Xu_Data{h,2};
    feature_num = 0;
    for i = 1:size(Amino_acids, 2)
        for j = 1:size(Amino_acids, 2)
            feature_num = feature_num + 1;
            s = strcat(Amino_acids{i}, '(?=....', Amino_acids{j},')'); % concatenate 2 amino acids. Dot is added to indicate there can be any amino acid in between (concatenate to come to this form for e.g. X(?=....X))
            match = regexp(sequence, s); % find the concatenated amino acids in the sequence/peptide
            Num = size(match, 2); % store the number of matches
            Xu_Data{h,feature_num + 1768} = Num; % save the number as its feature
            feature_label {1,feature_num + 1768} = s;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Second Phase: Calculate the Feature Rank using F-Score method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Number of +ve and -ve samples and Separating +ve and -ve samples

Num_of_Pos  = 0;
Num_of_Neg  = 0;
for i = 1:size(Xu_Data,1)
    if Xu_Data{i,3} == 1
       Num_of_Pos = Num_of_Pos + 1;
       Pos_sample(Num_of_Pos,:) = Xu_Data(i,:); 
    end
    if Xu_Data{i,3} == 0
       Num_of_Neg = Num_of_Neg + 1;
       Neg_sample(Num_of_Neg,:) = Xu_Data(i,:);
    end    
end
% Sum the columns of the Final_Data
for i = 5:size(Xu_Data,2)
   
    summ = 0;
    for j = 1:size(Xu_Data,1)
        
      summ = summ + Xu_Data{j,i};  
        
    end
    Feature_sum{1,i} = summ;
end

% Mean j-th feature in whole
for i = 5:size(Feature_sum, 2)

     Feature_mean_whole{1,i} = Feature_sum{1,i}/size(Xu_Data,1);    

end

% Mean j-th feature in positive

% Sum the columns of the Pos_sample
for i = 5:size(Pos_sample,2)
   
    summ = 0;
    for j = 1:size(Pos_sample,1)
        
      summ = summ + Pos_sample{j,i};  
        
    end
    Feature_sum_pos{1,i} = summ;
end
% Mean j-th feature in Pos_sample
for i = 5:size(Feature_sum_pos, 2)

     Feature_mean_pos{1,i} = Feature_sum_pos{1,i}/size(Pos_sample,1);    

end

% Mean j-th feature in negative

% Sum the columns of the Neg_sample
for i = 5:size(Neg_sample,2)
   
    summ = 0;
    for j = 1:size(Neg_sample,1)
        
      summ = summ + Neg_sample{j,i};  
        
    end
    Feature_sum_neg{1,i} = summ;
end
% Mean j-th feature in Neg_sample
for i = 5:size(Feature_sum_neg, 2)

     Feature_mean_neg{1,i} = Feature_sum_neg{1,i}/size(Neg_sample,1);    

end

% Denominator part of F_Score formula

% Pos part of deno
for i = 5:size(Pos_sample, 2)
    summ = 0;
    for j = 1:size(Pos_sample, 1)

        summ = summ + (Pos_sample{j,i} - Feature_mean_pos{1,i})^2;
        
    end
    Pos_deno{1,i} = (1/(size(Pos_sample, 1) - 1)) * summ;
end
% Neg part of deno
for i = 5:size(Neg_sample, 2)
    summ = 0;
    for j = 1:size(Neg_sample, 1)

        summ = summ + (Neg_sample{j,i} - Feature_mean_neg{1,i})^2;
        
    end
    Neg_deno{1,i} = (1/(size(Neg_sample, 1) - 1)) * summ;
end

% F-Score Calculation

for i = 5:size(Xu_Data, 2)
    
    if (Pos_deno{1,i} + Neg_deno{1,i}) == 0
        F_Score{1,i} = 0;
    
    else
        F_Score{1,i} = ((Feature_mean_pos{1,i}-Feature_mean_whole{1,i})^2 + (Feature_mean_neg{1,i}-Feature_mean_whole{1,i})^2)/(Pos_deno{1,i} + Neg_deno{1,i});
    end
end

Fscore_n_FeatureLabel = [F_Score; feature_label];

A = Fscore_n_FeatureLabel';
B = A(5:size(A,1),:);
C = sortrows(B,1);
Feature_Rank = flipud(C); % Feature ranks


