%%% Script to compare BN coverage within our work to Heide, H. et 2012 Cell
%%% Meta
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 


%Using INPARANOID: Eukaryotic Ortholog Groups, Release 8.0, December 2013
%Created by Nichollas Scott

% Import INPARANOID 
%comparing human to rat
fid_Human_to_rat=fopen('Output.H.sapiens-R.norvegicus.csv', 'rt');                %Input corum data base. It should have the gene names
                                                    % in the same format as
                                                    % the chromatograms in
                                                    % two columns
Import_Human_to_rat= textscan(fid_Human_to_rat,'%s\t', 'Delimiter',',');
No_human_to_rat=length(Import_Human_to_rat{1})/1;
Import_Human_to_rat = (reshape(Import_Human_to_rat{1,1},1,No_human_to_rat)');

%Extract rat to human entries
Human_to_rat_arry=cell(1,4);
counter=1;

%Human orthologs to rat proteins
for lookup=1:No_human_to_rat
       try %Error detector, ensure if uniprot string is found but not in the correct format it is ignored     
       orhtolog=deblank(Import_Human_to_rat(lookup));
       Split_ortholog=regexp(orhtolog,' ','split');
       Split_ortholog=[Split_ortholog{:,:}];
       Split_ortholog=strrep(Split_ortholog,' ','');
       Split_ortholog=Split_ortholog(~cellfun('isempty',Split_ortholog));
       Split_ortholog=strtrim(Split_ortholog);
       Split_ortholog_middle=regexp(Split_ortholog{2},'		','split');
       %divide up into array 
       Human_to_rat_arry(counter,1)=Split_ortholog(1); %Human accession number goes here
       Human_to_rat_arry(counter,2)=Split_ortholog_middle(1); %INPARANOID comparsion value here
       Human_to_rat_arry(counter,3)=Split_ortholog_middle(2); %Rat accession number goes here
       Human_to_rat_arry(counter,4)=Split_ortholog(3); %INPARANOID comparsion value here
       counter=counter+1; 
       end
end

%Replace empty cell with NaN
Human_to_rat_arry(cellfun(@isempty, Human_to_rat_arry))={'NaN'};



%% Import the data
[~, ~, raw] = xlsread('BNE_cell_meta_2012.xlsx','Sheet1');
raw(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),raw)) = {''};
cellVectors = raw(:,1);

%% Allocate imported array to column variable names
BNE_array = cellVectors(:,1);

%Divide up name to extract just the uniprot number
Split_BNE_array=regexp(BNE_array,' | ','split');
Number_entries_BNE=length(Split_BNE_array);

%Create Array to copy names to
BNE_unpacked_array=cell(Number_entries_BNE,1);

for BNE_counter=1:Number_entries_BNE
    Unpacked=Split_BNE_array{BNE_counter,1};
    
    %Copy uniprot name to array
    BNE_unpacked_array(BNE_counter)=Unpacked(1);
end

%mapped array
BNE_mapped_array=cell(Number_entries_BNE,2);

%Find human homologue of protein in BNE
for BNE_counter=1:Number_entries_BNE
    %look up protein in Inpariod array
    [col1 row1]=ind2sub(size(Human_to_rat_arry), strmatch(BNE_unpacked_array{BNE_counter},Human_to_rat_arry, 'exact'));
    
    if ~isempty(col1)
    BNE_mapped_array(BNE_counter,1)=Human_to_rat_arry(col1,1);
    BNE_mapped_array(BNE_counter,2)=Human_to_rat_arry(col1,3);
    else isempty(col1)
    BNE_mapped_array(BNE_counter,1)={'NaN'};
    BNE_mapped_array(BNE_counter,2)={'NaN'};
    end
end

%% Import the data
[~, ~, raw] = xlsread('LP_BNE_cell_meta_2012.xlsx','Sheet1');
raw(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),raw)) = {''};
cellVectors = raw(:,1);

%% Allocate imported array to column variable names
LP_BNE_array = cellVectors(:,1);

%Divide up name to extract just the uniprot number
Split_LP_BNE_array=regexp(LP_BNE_array,' | ','split');
Number_entries_LP_BNE=length(Split_LP_BNE_array);

%Create Array to copy names to
LP_BNE_unpacked_array=cell(Number_entries_LP_BNE,1);

for LP_BNE_counter=1:Number_entries_LP_BNE
    Unpacked=Split_LP_BNE_array{LP_BNE_counter,1};
    
    %Copy uniprot name to array
    LP_BNE_unpacked_array(LP_BNE_counter)=Unpacked(1);
end

%mapped array
LP_BNE_mapped_array=cell(Number_entries_LP_BNE,2);

%Find human homologue of protein in BNE
for LP_BNE_counter=1:Number_entries_LP_BNE
    %look up protein in Inpariod array
    [col1 row1]=ind2sub(size(Human_to_rat_arry), strmatch(LP_BNE_unpacked_array{LP_BNE_counter},Human_to_rat_arry, 'exact'));
    
    if ~isempty(col1)
    LP_BNE_mapped_array(LP_BNE_counter,1)=Human_to_rat_arry(col1,1);
    LP_BNE_mapped_array(LP_BNE_counter,2)=Human_to_rat_arry(col1,3);
    else isempty(col1)
    LP_BNE_mapped_array(LP_BNE_counter,1)={'NaN'};
    LP_BNE_mapped_array(LP_BNE_counter,2)={'NaN'};
    end
end



%% Import Majority Protein ID list
[~, Protein_IDs] = xlsread('Major Protein Groups.xlsx');


%Create an array showing overlap of proteins
BNE_Overlap_array=zeros(1,4);

%total number of entries in BNE
BNE_Overlap_array(1)=length(Protein_IDs);

%total number of entries in BNE
BNE_Overlap_array(2)=Number_entries_BNE;

%Find human homologue of protein in BNE
for BNE_counter=1:Number_entries_BNE
    
    %check if value is non Nan
    if ~isequal(BNE_mapped_array(BNE_counter,1),{'NaN'})
    
    [col1 row1]=ind2sub(size(Protein_IDs), strmatch(BNE_mapped_array(BNE_counter,1),Protein_IDs, 'exact'));
    if ~isempty(col1)
    %Donote shared
    BNE_Overlap_array(3)= BNE_Overlap_array(3)+1;  
    end
    else
    %Donote not shared
    BNE_Overlap_array(4)= BNE_Overlap_array(4)+1;   
    end
    
end


%Create an array showing overlap of proteins
LP_BNE_Overlap_array=zeros(1,4);

%total number of entries in BNE
LP_BNE_Overlap_array(1)=length(Protein_IDs);

%total number of entries in BNE
LP_BNE_Overlap_array(2)=Number_entries_LP_BNE;

%Find human homologue of protein in BNE
for LP_BNE_counter=1:Number_entries_LP_BNE
    
    %check if value is non Nan
    if ~isequal(LP_BNE_mapped_array(LP_BNE_counter,1),{'NaN'})
    
    [col1 row1]=ind2sub(size(Protein_IDs), strmatch(LP_BNE_mapped_array(LP_BNE_counter,1),Protein_IDs, 'exact'));
    if ~isempty(col1)
    %Donote shared
    LP_BNE_Overlap_array(3)= LP_BNE_Overlap_array(3)+1;  
    end
    else
    %Donote not shared
    LP_BNE_Overlap_array(4)= LP_BNE_Overlap_array(4)+1;   
    end
    
end

Output_name3=strcat('Summary_Venn_diagram_information.txt');

fid = fopen(Output_name3,'w');
fprintf (fid,'%s\t%s\t%s\t%s\t%s\t\n',...
'Proteins ID in BN-PCP-SILAC','Proteins in BN-PAGE Heide, H. et al 2012','Proteins in LP-BN-PAGE Heide, H. et al 2012','Overlap with  BN-PAGE Heide, H. et al 2012','Overlap with  LP-BN-PAGE Heide, H. et al 2012'); %Write Header
fprintf(fid,'%f	',BNE_Overlap_array(1));
fprintf(fid,'%f	',BNE_Overlap_array(2));
fprintf(fid,'%f	',LP_BNE_Overlap_array(2));
fprintf(fid,'%f	',BNE_Overlap_array(3));
fprintf(fid,'%f	',LP_BNE_Overlap_array(3));
fclose(fid);



%Create Venn diagram
f1=figure
subplot(2,2,1)
venn([BNE_Overlap_array(1) BNE_Overlap_array(2)],BNE_Overlap_array(3),'FaceColor',{myC(1,:),myC(2,:)},'FaceAlpha',{0.5,0.5},'EdgeColor','black')
title('Overlap of BN-PAGE');
legend('BN-PCP-SILAC', 'BN-PAGE, Heide, H. et al 2012','FontSize',6, 'Location', 'southoutside')
subplot(2,2,2)
venn([LP_BNE_Overlap_array(1) LP_BNE_Overlap_array(2)],LP_BNE_Overlap_array(3),'FaceColor',{myC(1,:),myC(2,:)},'FaceAlpha',{0.5,0.5},'EdgeColor','black')
title('Overlap of LP-BN-PAGE');
legend('BN-PCP-SILAC', 'BN-PAGE, Heide, H. et al 2012','FontSize',6, 'Location', 'southoutside')
Output_name5=strcat('Venn_diagram_BN_PAGE_overlap.pdf');
print(f1,'-dpdf','-r600', Output_name5);

