%Combined unique interaction list, this script is designed to read in Final
%interaction csv from the ROC script, generate a final combined unique
%interaction list and determine the global precision, recall, FDR of the
%dataset.

clear all

%Define Colours to use
myC= [30/255 144/255 255/255
      255/255 215/255 0/255
      178/255 34/255 34/255
      193/255 205/255 193/255];
  
colour_to_use=[0.254 0.411 0.882 %Colour: Royal Blue
              135/255 206/255 250/255 %Colour: Sky Blue
              0.28 0.82 0.8 %Colour: Medium Turquoise
              205/255 92/255 92/255 %Colour: Indian Red
              178/255 34/255 34/255 %Colour: Firebrick
              255/255 69/255 0 %Colour: Orange Red
              244/255 238/255 224/255 %Colour: Honeydew 2
              106/255 90/255 205/255 %Colour: Slate Blue
              34/255 139/255 34/255 %Colour: Forest Green
              222/255 184/255 135/255 %Colour: Burlywood
              186/255 85/255 211/255 %Colour: Medium Orchid
              219/255 112/255 147/255]; %Colour: Pale Violet Red 
          

%Import interaction list
Import_Interaction_70pc_nameA= 'Combined results_cyto\Final_Interactions_list_70_precision.csv';
fileID_A = fopen(Import_Interaction_70pc_nameA,'r');
Import_Interaction_70pc_1 = textscan(fileID_A, '%s%s%s%s%s%s%s%s%s%s%s%s%[^\n\r]', 'Delimiter', ',', 'HeaderLines' ,1, 'ReturnOnError', false);
fclose(fileID_A);

%Import mito list
Import_Interaction_70pc_nameB= 'Combined results_mito\Final_Interactions_list_70_precision.csv';
fileID_B = fopen(Import_Interaction_70pc_nameB,'r');
Import_Interaction_70pc_2 = textscan(fileID_B, '%s%s%s%s%s%s%s%s%s%s%s%s%[^\n\r]', 'Delimiter', ',', 'HeaderLines' ,1, 'ReturnOnError', false);
fclose(fileID_B);


%Count number of interactions
detected_interaction1=length(Import_Interaction_70pc_1{1,1});
detected_interaction2=length(Import_Interaction_70pc_2{1,1});

%Convert for cyto interaction to a list that can be used
interaction_A.unique_interactions=Import_Interaction_70pc_1{1,1};
interaction_A.proteinA =Import_Interaction_70pc_1{1,2};
interaction_A.proteinB =Import_Interaction_70pc_1{1,3};
interaction_A.centerA_formated_cyto=Import_Interaction_70pc_1{1,4};
interaction_A.centerB_formated_cyto=Import_Interaction_70pc_1{1,5};
interaction_A.Replicates_formated_cyto=Import_Interaction_70pc_1{1,6};
interaction_A.DeltaHeight_formated_cyto=Import_Interaction_70pc_1{1,7};
interaction_A.DeltaCenter_formated_cyto=Import_Interaction_70pc_1{1,8};
interaction_A.Deltawidth_formated_cyto=Import_Interaction_70pc_1{1,9};
interaction_A.DeltaEuc_formated_cyto=Import_Interaction_70pc_1{1,11};
interaction_A.proteinInCorum=Import_Interaction_70pc_1{1,11};
interaction_A.interactionInCorum=Import_Interaction_70pc_1{1,12};

%Convert for cyto interaction to a list that can be used
interaction_B.unique_interactions=Import_Interaction_70pc_2{1,1};
interaction_B.proteinA =Import_Interaction_70pc_2{1,2};
interaction_B.proteinB =Import_Interaction_70pc_2{1,3};
interaction_B.centerA_formated_mito=Import_Interaction_70pc_2{1,4};
interaction_B.centerB_formated_mito=Import_Interaction_70pc_2{1,5};
interaction_B.Replicates_formated_mito=Import_Interaction_70pc_2{1,6};
interaction_B.DeltaHeight_formated_mito=Import_Interaction_70pc_2{1,7};
interaction_B.DeltaCenter_formated_mito=Import_Interaction_70pc_2{1,8};
interaction_B.Deltawidth_formated_mito=Import_Interaction_70pc_2{1,9};
interaction_B.DeltaEuc_formated_mito=Import_Interaction_70pc_2{1,11};
interaction_B.proteinInCorum=Import_Interaction_70pc_2{1,11};
interaction_B.interactionInCorum=Import_Interaction_70pc_2{1,12};

%Determine the list of unique proteins for cyto and mito
unique_cyto=unique([interaction_A.proteinA,interaction_A.proteinB]);
unique_mito=unique([interaction_B.proteinA,interaction_B.proteinB]);


%Determine how many times each unique protein is observed to interact
%For cyto

%Create array to store numbers
Number_of_interaction=zeros(length(unique_cyto),1);

for unique_cyto_counter=1:length(unique_cyto)
   %reset number
   position_forward=0;
   position_reverse=0;
   
   %look up interaction cyto list
   [position_forward]=ind2sub(size(interaction_A.proteinA), strmatch(unique_cyto{unique_cyto_counter}, interaction_A.proteinA, 'exact')); 
   [position_reverse]=ind2sub(size(interaction_A.proteinB), strmatch(unique_cyto{unique_cyto_counter}, interaction_A.proteinB, 'exact'));     
   %Count how many interactions it was observed in
   Number_of_interaction(unique_cyto_counter)=(length(position_forward)+length(position_reverse));
  
end


%Create log log plot
%determine properties of distrubution
Histo_bins=[1,2,3,4,5,6,7,8,9,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,90];
his_data=histc(Number_of_interaction,round(linspace(1,max(Number_of_interaction),length(Histo_bins))));
Y_axis_cyto=log(his_data/length(unique_cyto));
X_axis_cyto=log(Histo_bins);
Figure1=figure
scatter(X_axis_cyto,Y_axis_cyto, 'MarkerEdgeColor','k',...
              'MarkerFaceColor',myC(2,:));
          
%Remove any value with inf
position_inf=isinf(Y_axis_cyto);
Y_axis_cyto(position_inf)=[];
X_axis_cyto(position_inf)=[];          

%calculate Polynomial curve that fits data
Curve_user_defined1=polyfit(X_axis_cyto',Y_axis_cyto,1);
polyfit_x_value1=round(Curve_user_defined1(1)*100)/100;
polyfit_consent1=round(Curve_user_defined1(2)*100)/100;
Equation_of_polyfit1=strcat('Equation Of The Line: ', mat2str(polyfit_x_value1),'x +',mat2str(polyfit_consent1));
%testing Correlation Coefficients
correlation_window=corrcoef(X_axis_cyto',Y_axis_cyto);
R_value1=round(correlation_window(2)*100000)/100000;
R_squared1=round((correlation_window(2)*correlation_window(2))*1000)/1000;
R_value_for_figure1=strcat('Correlation Coefficient, R : ', mat2str(R_value1));
R_squared_for_figure1=strcat('Coefficient of Determination, R^2 : ', mat2str(R_squared1));
 
k=Curve_user_defined1(1);
loga=Curve_user_defined1(2);
a=exp(loga);
hold on; 
plot(X_axis_cyto,k*X_axis_cyto+loga,':k')
hold on
plot(X_axis_cyto,k*X_axis_cyto+loga,':k')
hold on
plot(X_axis_cyto,k*X_axis_cyto+loga,':k')
legend('Data',sprintf('y=%.3f{}log(x)+log(%.3f)',k,a),...
sprintf(R_value_for_figure1),sprintf(R_squared_for_figure1)); 

      
%scatter(exp(linspace(1,0,20)),exp(linspace(0,1,20)))
%hold on
%plot()

xlim([0 max(log(Histo_bins))*1.1]);
ylim([-8 1.1]);
title('Vertex degree distribution of Cytoplasmic interactome', 'Fontsize',10);
xlabel('Vertex degrees','Fontsize',10);
ylabel('Fraction of the Proteins','Fontsize',10);
%Chnage X axis labels
set(gca, 'Xtick',[log(Histo_bins(:))],'XTickLabel',[Histo_bins],...
    'Ytick', [min(Y_axis_cyto) max(Y_axis_cyto) 1],....
    'YTickLabel',[round(min(his_data/length(unique_cyto))*10000)/10000 round(max(his_data/length(unique_cyto))*100)/100 1]);


print(Figure1,'-dpdf','-r600', 'Assessment_of_power_distrubatiuon_cytoplasm.pdf');


%Determine how many times each unique protein is observed to interact
%For Membrane
%Create array to store numbers
Number_of_interaction_membrane=zeros(length(unique_mito),1);

for unique_mito_counter=1:length(unique_mito)
   %reset number
   position_forward_membrane=0;
   position_reverse_membrane=0;
   
   %look up interaction cyto list
   [position_forward_membrane]=ind2sub(size(interaction_B.proteinA), strmatch(unique_mito{unique_mito_counter}, interaction_B.proteinA, 'exact')); 
   [position_reverse_membrane]=ind2sub(size(interaction_B.proteinB), strmatch(unique_mito{unique_mito_counter}, interaction_B.proteinB, 'exact'));     
   %Count how many interactions it was observed in
   Number_of_interaction_membrane(unique_mito_counter)=(length(position_forward_membrane)+length(position_reverse_membrane));
  
end

%Create log log plot
%determine properties of distrubution
Histo_bins=[1,2,3,4,5,6,7,8,9,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100,105,110];
his_data2=histc(Number_of_interaction_membrane,Histo_bins);
Y_axis_mito=log(his_data2/length(unique_mito));
X_axis_mito=log(Histo_bins);
Figure1=figure
scatter(X_axis_mito,Y_axis_mito, 'MarkerEdgeColor','k',...
              'MarkerFaceColor',myC(2,:));
          
%Remove any value with inf
position_inf=isinf(Y_axis_mito);
Y_axis_mito(position_inf)=[];
X_axis_mito(position_inf)=[];

%calculate Polynomial curve that fits data
Curve_user_defined1=polyfit(X_axis_mito',Y_axis_mito,1);
polyfit_x_value1=round(Curve_user_defined1(1)*100)/100;
polyfit_consent1=round(Curve_user_defined1(2)*100)/100;
Equation_of_polyfit1=strcat('Equation Of The Line: ', mat2str(polyfit_x_value1),'x +',mat2str(polyfit_consent1));
%testing Correlation Coefficients
correlation_window=corrcoef(X_axis_mito',Y_axis_mito);
R_value1=round(correlation_window(2)*100000)/100000;
R_squared1=round((correlation_window(2)*correlation_window(2))*1000)/1000;
R_value_for_figure1=strcat('Correlation Coefficient, R : ', mat2str(R_value1));
R_squared_for_figure1=strcat('Coefficient of Determination, R^2 : ', mat2str(R_squared1));
 
k=Curve_user_defined1(1);
loga=Curve_user_defined1(2);
a=exp(loga);
hold on; 
plot(X_axis_mito,k*X_axis_mito+loga,':k')
hold on
plot(X_axis_mito,k*X_axis_mito+loga,':k')
hold on
plot(X_axis_mito,k*X_axis_mito+loga,':k')
legend('Data',sprintf('y=%.3f{}log(x)+log(%.3f)',k,a),...
sprintf(R_value_for_figure1),sprintf(R_squared_for_figure1)); 

      
%scatter(exp(linspace(1,0,20)),exp(linspace(0,1,20)))
%hold on
%plot()

xlim([0 max(log(Histo_bins))*1.1]);
ylim([-8 1.1]);
title('Vertex degree distribution of membrane interactome', 'Fontsize',10);
xlabel('Vertex degrees','Fontsize',10);
ylabel('Fraction of the Proteins','Fontsize',10);
%Chnage X axis labels
set(gca, 'Xtick',[log(Histo_bins(:))],'XTickLabel',[Histo_bins],...
    'Ytick', [min(Y_axis_mito) max(Y_axis_mito) 1],....
    'YTickLabel',[round(min(his_data2/length(unique_mito))*10000)/10000 round(max(his_data2/length(unique_mito))*100)/100 1]);


print(Figure1,'-dpdf','-r600', 'Assessment_of_power_distrubatiuon_membrane.pdf');





