%%% Using GOCC count the number of membrane interactions
%Created by Nichollas Scott
clear all

%% colour plattes to outputs
%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 
              185/255 211/255 238/255
              138/255 43/255 226/255  
              238/255 180/255 180/255
              205/255 51/255 51/255
              238/255 118/255 33/255
              143/255 188/255 143/255
              84/255 255/255 159/255
              238/255 230/255 133/255]; 

%Import Annotation GOCC assignments from Persaus, remove Proteingroups
%without assignment for simplicity
fid_annotation=fopen('Annotation_protein_group_information_GOCC.csv', 'rt');                
                                                   
Import_annotation= textscan(fid_annotation,'%s\t','Delimiter','\n');
No_of_cells=length(Import_annotation{1});

%Remove " from files
Import_annotation=strrep(Import_annotation{:},'"','');

%create array of GOCC
for counter=1:No_of_cells
    seperated_values=regexp(Import_annotation(counter,1),'(;|,)','split');
    Split_values=[seperated_values{:,:}];
    %determine length of array
    row_number=length(Split_values);
    %Copy value to array
    for counter2=1:row_number
    Imported_seperated(counter,counter2)=Split_values(counter2);
    end
end

%replace empty cell with NaN
emptyCells = cellfun(@isempty,Imported_seperated);
Imported_seperated(emptyCells) = {'NaN'};


%Import interaction list
Import_Interaction_70pc= importdata('Final_Interactions_list_70_precision.csv');

%Count number of interactions
detected_interaction=length(Import_Interaction_70pc.textdata);

%split interaction, note removed title
for counter=2:detected_interaction
    seperated_values=regexp(Import_Interaction_70pc.textdata(counter,1),'_','split');
    Split_values=[seperated_values{:,:}];

    %determine length of array
    row_number=length(Split_values);
    %Copy value to array
    for counter2=1:row_number
    Imported_seperated2(counter-1,counter2)=Split_values(counter2);
    end

end

%Look up proteinID for both interactions and test if 
% A) the proteinID have GOCC
% B) the GOCC are referred to as membrane
% C) the GOCC are the same between the proteinIDs

detected_interaction=length(Imported_seperated2);

%Create array 
contains_membrane=zeros(detected_interaction,2);
contains_membrane_terms=cell(detected_interaction,1);
Terms_A=cell(1,1);
Terms_B=cell(1,1);

for counter=1:detected_interaction

%look up if protein is part of human corum database
[position1]=ind2sub(size(Imported_seperated), strmatch(Imported_seperated2{counter,1}, Imported_seperated, 'exact'));

%look up if protein is part of human corum database
[position2]=ind2sub(size(Imported_seperated), strmatch(Imported_seperated2{counter,2}, Imported_seperated, 'exact'));

%Find if GOCC contains the word "membrane"
GOCC_tocheck1=[Imported_seperated(position1,:)];
%Save data to use to extract terms
GOCC_tocheck1B=GOCC_tocheck1;
GOCC_tocheck1=strfind(GOCC_tocheck1,'membrane');
emptyCells1 = ~cellfun(@isempty,GOCC_tocheck1);
index1=find(emptyCells1>0);

%Save list of GOCC terms for analysis
interal_counterA=1;
for term_counter=2:length(GOCC_tocheck1)
    if ~isequal(GOCC_tocheck1B(term_counter),{'NaN'})
     Terms_A(counter,interal_counterA)=GOCC_tocheck1B(term_counter);
     interal_counterA=1+interal_counterA;
    end
end

%Find if GOCC contains the word "membrane"
GOCC_tocheck2=[Imported_seperated(position2,:)];
%Save data to use to extract terms
GOCC_tocheck2B=GOCC_tocheck2;
GOCC_tocheck2=strfind(GOCC_tocheck2,'membrane');
emptyCells2 = ~cellfun(@isempty,GOCC_tocheck2);
index2=find(emptyCells2>0);

%Save list of GOCC terms for analysis
interal_counterB=1;
for term_counter=2:length(GOCC_tocheck2)
    if ~isequal(GOCC_tocheck2B(term_counter),{'NaN'})
     Terms_B(counter,interal_counterB)=GOCC_tocheck2B(term_counter);
     interal_counterB=1+interal_counterB;
    end
end

if ~isempty(index1) | ~isempty(index2)
    %Write out a number
    contains_membrane(counter,1)=1;
end

 %compare GOCC terms and record which ones are the same
 ProteinA= Imported_seperated(position1,index1);
 ProteinB= Imported_seperated(position2,index2);
 
 %Determine which terms are shared
 GOCC_terms_shared=ismember(ProteinA, ProteinB);
 
 %Count how many terms are shared
 contains_membrane(counter,2)=sum(GOCC_terms_shared);
 
 %write out terms
 contains_membrane_terms(counter,1)={strjoin(ProteinA(GOCC_terms_shared),' ; ')};

end

%Write file showing list of homologs
fid = fopen('Membrane_interactions.csv','w');
fprintf (fid,'%s,%s,%s,%s,%s,%s,%s,%s,%s\n',... %header for output
'Interaction','Protein A','Protein B',...
'Center A','Center B','Replicate','GOOC term containing the word "Membrane"',...
'Number of GOCC terms shared','GOCC terms shared by interactors'); %Write Header
for write_out_loop2=1:detected_interaction
fprintf(fid,'%s,',Import_Interaction_70pc.textdata{write_out_loop2+1,1});
fprintf(fid,'%s,',Import_Interaction_70pc.textdata{write_out_loop2+1,2});
fprintf(fid,'%s,',Import_Interaction_70pc.textdata{write_out_loop2+1,3});
fprintf(fid,'%s,',Import_Interaction_70pc.textdata{write_out_loop2+1,4});
fprintf(fid,'%s,',Import_Interaction_70pc.textdata{write_out_loop2+1,5});
fprintf(fid,'%s,',Import_Interaction_70pc.textdata{write_out_loop2+1,6});
fprintf(fid,'%f,',contains_membrane(write_out_loop2,1));
fprintf(fid,'%f,',contains_membrane(write_out_loop2,2));
fprintf(fid,'%s,',contains_membrane_terms{write_out_loop2,1});
fprintf(fid,',\n');
end
fclose(fid);


%Visualization of interactions based on terms

%pie chart of number of interactions with exact same Go term, related go
%term (contains the work membrane in both), doesnt contain the work
%membrane in one

Stats(1,1)=sum(contains_membrane(:,1)~=0);
Stats(1,2)=sum(contains_membrane(:,2)~=0);
Stats(1,3)=numel(contains_membrane(:,1));


stats_pc(1,1)=Stats(1,2)/Stats(1,3);
stats_pc(1,2)=(Stats(1,1)-Stats(1,2))/Stats(1,3);
stats_pc(1,3)=1-(stats_pc(1,1)+stats_pc(1,2));

%Membrane Go terms
f1 = figure;
h=pie(stats_pc)
hp = findobj(h, 'Type', 'patch');
set(hp(1), 'FaceColor', myC(2,:));
set(hp(2), 'FaceColor', colour_to_use(2,:));
set(hp(3), 'FaceColor', colour_to_use(4,:));
legend('GOCC term containing the word "Membrane"','GOCC terms shared by interactors','GOCC term not shared',6, 'Location', 'southoutside')
title('Membrane GOCC terms');
print(f1,'-dpdf','-r600', 'GOCC membrane pie graph');
print(f1,'-dpng','-r600', 'GOCC membrane pie graph.png');

%Combined terms
Combined_terms=[Terms_A,Terms_B];
%remove empty cells
Combined_terms_new = Combined_terms(~cellfun(@isempty, Combined_terms));
unique_terms=unique(Combined_terms_new);
Num_unique_terms=length(unique_terms);

%Count terms
Matrix_terms=zeros(Num_unique_terms,Num_unique_terms);
Matrix_terms_specific=zeros(Num_unique_terms,Num_unique_terms);

%length term A and B
Length_termA=size(Terms_A);
Length_termB=size(Terms_B);


for term_counter=1:Num_unique_terms
%Write out text to keep track of which term is being processed    
strcat('now processing term_',mat2str(term_counter),'_of_',mat2str(Num_unique_terms))

for counterA=2:detected_interaction
for counterB=1:Length_termA(2)

%determine terms is observed in protein A

if ~cellfun('isempty',unique_terms(term_counter))
if isequal(Terms_A(counterA,counterB),unique_terms(term_counter))
 

%term is observed count how many times each of the terms were see in B
for counter2B=1:Length_termB(2)
%create varible to avoid overhead
test_varible=Terms_B(counterA,counter2B);

%Determine if variable is observed
for term_counterB=1:Num_unique_terms

if ~cellfun('isempty',unique_terms(term_counterB))
if isequal(test_varible,unique_terms(term_counterB))

%Count how many terms co appear    
Matrix_terms(term_counter,term_counterB)=Matrix_terms(term_counter,term_counterB)+1;
Matrix_terms_specific(term_counter,term_counterB)=1;    
end

end
end
end

end
end

end
end

end


%Save Matrix file to be used
%Save as .mat to load quickly
save('Matrix_GOCC_terms.mat','Matrix_terms');
save('Matrix_terms_specific.mat','Matrix_terms_specific');

%load mat
%load('Matrix_GOCC_terms.mat','Matrix_terms')

%Save csv file 
csvwrite('Matrix_GOCC_terms.csv',Matrix_terms);
csvwrite('Matrix_Specific_GOCC_terms.csv',Matrix_terms_specific);



%Sum Matrix and generate a histogram
max_matrix_terms=max(Matrix_terms);
max_matrix_terms_specific=max(Matrix_terms_specific);

sum_matrix_terms=sum(Matrix_terms);
sum_matrix_terms_specific=sum(Matrix_terms_specific);

%sort
[sorted_max_matrix_terms, sort_order1]=sort(max_matrix_terms);
[sorted_sum_matrix_terms, sort_order2]=sort(sum_matrix_terms);
[sorted_max_matrix_terms_specific, sort_order3]=sort(max_matrix_terms_specific);
[sorted_sum_matrix_terms_specific, sort_order4]=sort(sum_matrix_terms_specific);


%Re-order matrix
re_ordered_matrix=Matrix_terms(:,sort_order1);
re_ordered_matrix2=Matrix_terms_specific(:,sort_order4);

%Show figure of number of unique GOCC terms
f2=figure;
%plot 1
subplot(1,5,1)
bar(re_ordered_matrix', 'stacked');
axis([0 length(sorted_sum_matrix_terms) 0 max(sorted_sum_matrix_terms)*1.05]);
title('All GOCC terms');
ax = gca;
ax.XTickLabelRotation = -90;
h = get(gca, 'xlabel');
set(h, 'FontSize', 8);

%plot 2
subplot(1,5,2:5)
hb = bar(re_ordered_matrix(:,(length(sorted_sum_matrix_terms)-29):length(sorted_sum_matrix_terms))', 'stacked');
axis([0 31 0 max(sorted_sum_matrix_terms)*1.05]);
title('Top 30 terms');
ax = gca;
ax.XTick=[1:30];
ax.XTickLabel=[(length(sorted_sum_matrix_terms)-29):length(sorted_sum_matrix_terms)];
ax.XTickLabelRotation = -90;
h = get(gca, 'xlabel');
set(h, 'FontSize', 8);
set(gca,'ytick',[]);
set(gca,'yticklabel',[]);

% set super title
suptitle('All GOCC terms combinations observed within interactions');
print(f2,'-dpdf','-r600', 'GOCC_histogram_terms_every_combination.pdf');
%Save as a non-vector graphic
print(f2,'-dpng','-r600', 'GOCC_histogram_terms_every_combination.png');

%plot 3, figure legend
f4A=figure;
hb = bar(re_ordered_matrix(:,(length(sorted_sum_matrix_terms)-29):length(sorted_sum_matrix_terms))', 'stacked');
hl=legend(hb((length(sorted_sum_matrix_terms)-29):length(sorted_sum_matrix_terms)),...
unique_terms{sort_order1((length(sorted_sum_matrix_terms)-29):length(sorted_sum_matrix_terms))},'Location','EastOutside','FontSize',6);
set(hb, 'visible', 'off');
set(gca, 'visible', 'off');
print(f4A,'-dpdf','-r600', 'legend_GOCC_histogram_terms_every_combinations.pdf')
close(f4A)


%Write out summary text file
fid = fopen('All GOCC terms combinations observed within interactions.txt','w');
fprintf (fid,'%s\t%s\t\n',... %header for output
'GOCC term','Number of observed interations'); %Write Header
for write_out_loop2=1:length(sorted_sum_matrix_terms)
%Write out data            
fprintf(fid,'%s	',unique_terms{sort_order2(write_out_loop2)});
fprintf(fid,'%f	',sorted_sum_matrix_terms(write_out_loop2));
fprintf(fid,'\n');
end
fclose(fid);

%Show figure of number of unique GOCC terms
f3=figure;
%plot 1
subplot(1,5,1)
bar(sorted_sum_matrix_terms_specific);
axis([0 length(sum_matrix_terms_specific) 0 max(sum_matrix_terms_specific)*1.05]);
title('All GOCC terms','FontSize', 8);
ax = gca;
ax.XTickLabelRotation = -90;
h = get(gca, 'xlabel');
set(h, 'FontSize', 8);

%plot 2
subplot(1,5,2:5)
hb = bar(sorted_sum_matrix_terms_specific((length(sum_matrix_terms_specific)-29):length(sum_matrix_terms_specific)));
axis([0 31 0 max(sum_matrix_terms_specific)*1.05]);
title('Top 30 terms','FontSize', 8);
ax = gca;
ax.XTick=[1:30];
ax.XTickLabel=[(length(sum_matrix_terms_specific)-29):length(sum_matrix_terms_specific)];
ax.XTickLabelRotation = -90;
h = get(gca, 'xlabel');
set(h, 'FontSize', 8);
set(gca,'ytick',[]);
set(gca,'yticklabel',[]);

% set super title
suptitle('Unique GOCC terms combinations observed within interactions');
print(f3,'-dpdf','-r600', 'GOCC_histogram_terms_unique_combinations.pdf')
%Save as a non-vector graphic
print(f3,'-dpng','-r600', 'GOCC_histogram_terms_unique_combination.png');
close(f3)

%plot 3, figure legend
f4B=figure;
hb = bar(re_ordered_matrix2(:,(length(sum_matrix_terms_specific)-19):length(sum_matrix_terms_specific))');
hl=legend(hb((length(sum_matrix_terms_specific)-19):length(sum_matrix_terms_specific)),...
{unique_terms{sort_order1((length(sum_matrix_terms_specific)-19):length(sum_matrix_terms_specific))}},'FontSize',6);
set(hb, 'visible', 'off');
set(gca, 'visible', 'off');
print(f4B,'-dpdf','-r600', 'legend_GOCC_histogram_terms_unique_combinations.pdf')
close(f4B)

%Write out summary text file
fid = fopen('Unique GOCC terms combinations observed within interactions.txt','w');
fprintf (fid,'%s\t%s\t\n',... %header for output
'GOCC term','Number of observed interations'); %Write Header
for write_out_loop2=1:length(sorted_sum_matrix_terms)
%Write out data            
fprintf(fid,'%s	',unique_terms{sort_order4(write_out_loop2)});
fprintf(fid,'%f	',sorted_sum_matrix_terms_specific(write_out_loop2));
fprintf(fid,'\n');
end
fclose(fid);

%Write out summary text file of GOCC stats
Average_number_of_GOCC_terms_combinations=mean(sorted_sum_matrix_terms_specific);
stdev_of_GOCC_terms_combinations=std(sorted_sum_matrix_terms_specific);
median_of_GOCC_terms_combinations=median(sorted_sum_matrix_terms_specific);


%Write out summary text file
fid = fopen('Summary of GOCC terms stats.txt','w');
fprintf (fid,'%s\t%s\t%s\t\n',... %header for output
'Average number of GOCC terms','std dev of GOCC terms','mean number of GOCC terms'); %Write Header
%Write out data            
fprintf(fid,'%f	',Average_number_of_GOCC_terms_combinations(1));
fprintf(fid,'%f	',stdev_of_GOCC_terms_combinations(1));
fprintf(fid,'%f	',median_of_GOCC_terms_combinations(1));
fclose(fid);

