%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%              ROC analysis of PCP-SEC-SILAC Datasets
%%   Created by Anders Kristensen, modified by Nichollas Scott, 
%%                             August 2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ensure: strjoin.m, myPCA.m, myCenter.m, myWhiten.m, d2p.m and tsne.m are 
%located in the parent dictory

%% ROC curves SEC-PCP-SILAC
clear all
matlabpool close force local
matlabpool;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% global varibles 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%User defined setting
number_of_replicates=3;
number_of_channels=2;

%Add scripts within this dictory to functions that can be called
addpath(pwd);

%% Import Majority Protein ID list
[~, Protein_IDs] = xlsread('Major_protein_groups.xlsx');
%Maxqaunt identified Majority protein ID for each Protein Group

%% Import Corum interactions 
fid=fopen('Mapped_human_Corum_list_20150109.csv', 'rt');    %Input corum data base.
                                                                % in the same format as
                                                                % the chromatograms in
                                                                % two columns
Corum_Import= textscan(fid,'%s\t', 'Delimiter',',');
fclose(fid);
No=length(Corum_Import{1})/2;
Corum_Protein_names = (reshape(Corum_Import{1,1},2,No)');
Unique_Corum = unique(Corum_Protein_names);

%% Import Corum Complex list, 
[~, ~, Corum_complexes] = xlsread('Corum_2012_human.xlsx','Sheet1','A2:L1847');
Corum_complexes(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),Corum_complexes)) = {''};
cellVectors = Corum_complexes(:,[2,3,4,5,7,9,10,11,12]);
Corum_complexes = Corum_complexes(:,[1,6,8]);

%% Replace non-numeric cells with NaN
R = cellfun(@(x) ~isnumeric(x) && ~islogical(x),Corum_complexes); % Find non-numeric cells
Corum_complexes(R) = {NaN}; % Replace non-numeric cells

%% Create output variable
Corum_complexes_data = reshape([Corum_complexes{:}],size(Corum_complexes));

%% Import replicate data 
%define Raw SILAC ratios data, This is the output from the alignment script
MvsL_filename_Raw_rep1='MvsL_alignment\Realignment\Adjusted_MvsL_Raw_for_ROC_analysis_rep1.csv';
MvsL_filename_Raw_rep2='MvsL_alignment\Realignment\Adjusted_MvsL_Raw_for_ROC_analysis_rep2.csv';
MvsL_filename_Raw_rep3='MvsL_alignment\Realignment\Adjusted_MvsL_Raw_for_ROC_analysis_rep3.csv';
HvsL_filename_Raw_rep1='HvsL_alignment\Realignment\Adjusted_HvsL_Raw_for_ROC_analysis_rep1.csv';
HvsL_filename_Raw_rep2='HvsL_alignment\Realignment\Adjusted_HvsL_Raw_for_ROC_analysis_rep2.csv';
HvsL_filename_Raw_rep3='HvsL_alignment\Realignment\Adjusted_HvsL_Raw_for_ROC_analysis_rep3.csv';

MvsL_filename_gaus_rep1='MvsL_alignment\Realignment\Adjusted_Combined_OutputGaus_rep1.csv';
MvsL_filename_gaus_rep2='MvsL_alignment\Realignment\Adjusted_Combined_OutputGaus_rep2.csv';
MvsL_filename_gaus_rep3='MvsL_alignment\Realignment\Adjusted_Combined_OutputGaus_rep3.csv';
HvsL_filename_gaus_rep1='HvsL_alignment\Realignment\Adjusted_Combined_OutputGaus_rep1.csv';
HvsL_filename_gaus_rep2='HvsL_alignment\Realignment\Adjusted_Combined_OutputGaus_rep2.csv';
HvsL_filename_gaus_rep3='HvsL_alignment\Realignment\Adjusted_Combined_OutputGaus_rep3.csv';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%create list of file names
%MvsL
List_of_Raw_filename={MvsL_filename_Raw_rep1,MvsL_filename_Raw_rep2,MvsL_filename_Raw_rep3,...
                      HvsL_filename_Raw_rep1,HvsL_filename_Raw_rep2,HvsL_filename_Raw_rep3};
List_of_Gaus_filename={MvsL_filename_gaus_rep1,MvsL_filename_gaus_rep2,MvsL_filename_gaus_rep3,...
                      HvsL_filename_gaus_rep1,HvsL_filename_gaus_rep2,HvsL_filename_gaus_rep3};


%% Allocate imported array to column variable names
Corum.Complexid = Corum_complexes_data(:,1);
Corum.Complexname = cellVectors(:,1);
Corum.Synonyms = cellVectors(:,2);
Corum.organism = cellVectors(:,3);
Corum.subunitsUniProtIDs = cellVectors(:,4);
Corum.subunitsEntrezIDs = Corum_complexes_data(:,2);
Corum.proteincomplexpurificationmethod = cellVectors(:,5);
Corum.PubMedid = Corum_complexes_data(:,3);
Corum.FunCatcategories = cellVectors(:,6);
Corum.functionalcomment = cellVectors(:,7);
Corum.diseasecomment = cellVectors(:,8);
Corum.subunitcomment = cellVectors(:,9);

%Split UniprotID to individual entries
length_of_corum=length(Corum.subunitsUniProtIDs);
Corum.subunitsUniProtIDs_split=cell(length_of_corum,1);
Split_Complex=regexp(Corum.subunitsUniProtIDs,',','split');
for Complex_no=1:length_of_corum
%using reg expression split entries in corum
number_proteins_in_complex=length(Split_Complex{Complex_no});
if number_proteins_in_complex > 1
Complex_to_write_out=Split_Complex{Complex_no};
    for write_out_counter1=1:number_proteins_in_complex
        Corum.subunitsUniProtIDs_split{Complex_no,write_out_counter1}=Complex_to_write_out(1,write_out_counter1); 
    end
elseif number_proteins_in_complex == 1
    Corum.subunitsUniProtIDs_split{Complex_no,1}=Split_Complex{Complex_no};
end
end

%Create file to save replicate information to
Save_replicate_folder_name=cell((number_of_replicates*number_of_channels),1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Replicate processing loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for replicate_counter=1:(number_of_replicates*number_of_channels)
    
% write out message to user to tell them what replicate is being processed 
strcat('now processing replicate ',mat2str(replicate_counter)) 
replicate_folder_name=strcat('Replicate ',mat2str(replicate_counter));
Save_replicate_folder_name{replicate_counter}=replicate_folder_name;


mkdir(Save_replicate_folder_name{replicate_counter});

%% Import MvsL
%Raw data for analysis of euclidean distance
%import data files summary data
Maxquant = fopen (List_of_Raw_filename{replicate_counter});
Maxquant_raw = textscan(Maxquant, '%s', 'Delimiter',',');
fclose(Maxquant);

%Reshape data for use in analysis
Raw_Dimension=size(Maxquant_raw{:});
Raw_data_reshaped=reshape(Maxquant_raw{:},68,(Raw_Dimension(1)/68))';
Raw_data_reshaped_size=size(Raw_data_reshaped);

%divide up reshaped data in correct form
num_val=Raw_data_reshaped(2:end,3:end);
num_val=cellfun(@str2num,num_val);

%Gaus data data for analysis of center analysis
Gaus_import_name = fopen (List_of_Gaus_filename{replicate_counter});
Gaus_import = textscan(Gaus_import_name, '%s', 'Delimiter',',');
fclose(Gaus_import_name);

%Reshape data for use in analysis
Gaus_import_Dimension=size(Gaus_import{:});
Gaus_import_reshaped=reshape(Gaus_import{:},7,(Gaus_import_Dimension(1)/7))';
Gaus_import_reshaped_size=size(Gaus_import_reshaped);

%chnage director to write replicate information into its own folder
cd(replicate_folder_name);

%create directories
mkdir('ROC analysis related Figures');
mkdir('ROC analysis matrixes');
mkdir('ROC analysis data');

%% Pre-process data to remove Gaussian above fraction five
%Create array to use for comparsions
H_raw =Gaus_import_reshaped(2:end,2);
C_raw =Gaus_import_reshaped(2:end,3);
W_raw =Gaus_import_reshaped(2:end,4);
%Covert to number
H_raw=cellfun(@str2num,H_raw);
C_raw=cellfun(@str2num,C_raw);
W_raw=cellfun(@str2num,W_raw);

Len_raw=size(H_raw);

%Filter data to remove NaN values
num_val(isnan(num_val))= 0.05;

%% The Chromatograms
Chromatograms_raw = num_val;
%reset Gaussian_to_remove
Gaussian_to_remove=zeros(0,0);

%find all guassian with center within the first five fraction
Gaussian_to_remove_counter=1;
for Filter_fractions=1:Len_raw(1)
    if C_raw(Filter_fractions)<5
        Gaussian_to_remove(Gaussian_to_remove_counter)=Filter_fractions;
        Gaussian_to_remove_counter=Gaussian_to_remove_counter+1;
    end
end

%remove gaussian with centers below five
H_raw([Gaussian_to_remove(:)])=[];
C_raw([Gaussian_to_remove(:)])=[];
W_raw([Gaussian_to_remove(:)])=[];
Chromatograms_raw(Gaussian_to_remove(:),:)=[];

H=H_raw;
C=C_raw;
W=W_raw;
Chromatograms=Chromatograms_raw;
%Calculate size of matrix
Len=size(H);
Chromatograms_dimension=size(Chromatograms);

%% Generate array of points to assess for normalised Gaussian
points_to_asses=0.25:0.25:100;
gaussian_fit=zeros(Len(1),length(points_to_asses));

%Generate normalised Gaussian
for row_counter=1:Len(1)
for test_value=1:length(points_to_asses)    

%generate normalised Gaus profile
gaussian_fit(row_counter,test_value)=(1*exp(-((points_to_asses(test_value)- C_raw(row_counter))...
/(W_raw(row_counter))).^2));
end
end

%clean data to remove vaues below  <0.1
gaussian_fit(gaussian_fit<0.001)=0;


%Calculate distances
Dist.Euc = squareform(pdist(Chromatograms,'euclidean'));
Dist.Center = squareform(pdist(C,'euclidean'));
Dist.Height = squareform(pdist(H,'euclidean'));
Dist.Width = squareform(pdist(W,'euclidean'));
Dist.Gaussian_fits= squareform(pdist(gaussian_fit,'euclidean'));

%% Combine Majority Protein ID with Proteins Identified in replicate
Protein.Isoform=Gaus_import_reshaped(2:end,1);
Protein.Height=cellfun(@str2num,Gaus_import_reshaped(2:end,2));
Protein.Width=cellfun(@str2num,Gaus_import_reshaped(2:end,4));
Protein.Center=cellfun(@str2num,Gaus_import_reshaped(2:end,3));

%remove protein name with center less then 5
Protein.Isoform(Gaussian_to_remove(:))=[];
Protein.Height(Gaussian_to_remove(:))=[];
Protein.Width(Gaussian_to_remove(:))=[];
Protein.Center(Gaussian_to_remove(:))=[];

%Create unique protein ID with center and protein name
for unique_counter=1:Len(1)
Protein.Unique_ID{unique_counter,1}=strcat(mat2str(Protein.Center(unique_counter)),'_',Protein.Isoform{unique_counter});
end

Dimensions_Gaus_import= length(Protein.Isoform);
Dimension_of_Protein_IDs=size(Protein_IDs);

%Convert the Protein information with isoform information to Structure with
%both isoform and no isoform information
for Gaus_import_counter1 = 1:Dimensions_Gaus_import(1)
    Find_Gaus_isoforms=strfind(Protein.Isoform(Gaus_import_counter1), '-');
    if ~isempty(Find_Gaus_isoforms{1}) == 1
        length_of_Find_Gaus_Protein=length(Protein.Isoform{Gaus_import_counter1});
        number_of_characters_to_remove2=Find_Gaus_isoforms{1}-1;
        Protein.NoIsoform(Gaus_import_counter1,1)={Protein.Isoform{Gaus_import_counter1}(1:number_of_characters_to_remove2)};
    elseif isempty(Find_Gaus_isoforms{1}) == 1
        Protein.NoIsoform(Gaus_import_counter1,1)=Protein.Isoform(Gaus_import_counter1);
    end
end

%Process Majority Protein ID data
for Protein_ID_counter1 =2:Dimension_of_Protein_IDs(1)
    for Protein_ID_counter2 =1:Dimension_of_Protein_IDs(2)-sum(cellfun('isempty',Protein_IDs(Protein_ID_counter1,:)))
    Find_isoforms=strfind(Protein_IDs(Protein_ID_counter1,Protein_ID_counter2), '-');
    if ~isempty(Find_isoforms{1}) == 1
        length_of_Protein_ID=length(Protein_IDs{Protein_ID_counter1,Protein_ID_counter2});
        number_of_characters_to_remove=Find_isoforms{1}-1;
        Protein_IDS_no_isoform(Protein_ID_counter1,Protein_ID_counter2)=cellstr(Protein_IDs{Protein_ID_counter1,Protein_ID_counter2}(1:number_of_characters_to_remove));
    elseif isempty(Find_isoforms{1}) == 1
        Protein_IDS_no_isoform(Protein_ID_counter1,Protein_ID_counter2)=Protein_IDs(Protein_ID_counter1,Protein_ID_counter2);
    end
    end
end

%Remove non unique entries on each row for MajorIDs
for Protein_ID_counter3 =2:Dimension_of_Protein_IDs(1)
     Array_to_check1=Protein_IDS_no_isoform(Protein_ID_counter3,:);
     Array_to_check_empty=cellfun(@isempty,Array_to_check1);
     Array_to_check1(Array_to_check_empty)=[]; % Replace empty cells with NaN
     uniqueCells=unique(Array_to_check1);
     parfor Protein_ID_counter4 =1:length(uniqueCells) 
            Protein_IDS_no_isoform_no_dup_processed{Protein_ID_counter3,Protein_ID_counter4}=uniqueCells(Protein_ID_counter4);
     end
end

%Define the dimensions of MajorIDs with duplicates removed
no_dup_processed_dimensions=size(Protein_IDS_no_isoform_no_dup_processed);

%Create strcuature array to fill
Protein.MajorIDs=cell(Dimensions_Gaus_import,Dimension_of_Protein_IDs(2));  
Protein.MajorID_NoIsoforms=cell(Dimensions_Gaus_import,Dimension_of_Protein_IDs(2));

%Look up Major Protein group Information
% this loop looks up the position of the Protein observed within the replicate
% within the MajorID list MajorIDs with isoforms and populates the Protein
% list with all the isoform information from the Major Protein ID list
for Gaus_import_counter1 = 1:Dimensions_Gaus_import(1)
    Lookup_protein_name=Protein.Isoform(Gaus_import_counter1);
    for  Check_counter1 = 1:Dimension_of_Protein_IDs(1)
        for  Check_counter2 = 1:Dimension_of_Protein_IDs(2)
           if isequal(Lookup_protein_name(1),Protein_IDs(Check_counter1,Check_counter2))
               %define what row to look up
               Row_within_major_ID(Gaus_import_counter1,1)=Check_counter1(1);
               
               %copy names to MajorID from Protein_IDs raw data
               Array_to_check1=Protein_IDs(Row_within_major_ID(Gaus_import_counter1,1),:);
               Array_to_check_empty1=cellfun(@isempty,Array_to_check1);
               
               for Check_counter3 = 1:nnz(~Array_to_check_empty1)
               Protein.MajorIDs{Gaus_import_counter1,Check_counter3}={Protein_IDs{Row_within_major_ID(Gaus_import_counter1,1),Check_counter3}};
               end
               
               %copy names to MajorID_NoIsoforms from Protein_IDS_no_isoform_no_dup_processed
               Array_to_check2=Protein_IDS_no_isoform_no_dup_processed(Row_within_major_ID(Gaus_import_counter1,1),:);
               Array_to_check_empty2=cellfun(@isempty,Array_to_check2);
               
               for Check_counter4 = 1:nnz(~Array_to_check_empty2)
               Protein.MajorID_NoIsoforms{Gaus_import_counter1,Check_counter4}=Protein_IDS_no_isoform_no_dup_processed{Row_within_major_ID(Gaus_import_counter1,1),Check_counter4};
               end
           end
        end
    end
end

%Create a 1D array containing all Protein.MajorIDs and Protein.MajorIDs_Noisoform with row numbers
Protein_location_array.MajorID_NoIsoforms=cell(1,1);
Protein_location_array.MajorID=cell(1,1);
Protein_location_array.row=cell(1,1);
location_counter_1= 1;
for Isoform_counter1= 1:Dimensions_Gaus_import(1)
    for Isoform_counter2= 1:sum(~cellfun('isempty',Protein.MajorID_NoIsoforms(Isoform_counter1,:)))
        Protein_location_array.MajorID_NoIsoforms(location_counter_1,1)=Protein.MajorID_NoIsoforms(Isoform_counter1,Isoform_counter2);
        Protein_location_array.MajorID(location_counter_1,1)=Protein.MajorIDs(Isoform_counter1,Isoform_counter2);
        Protein_location_array.row{location_counter_1,1}=Isoform_counter1;
        location_counter_1=1+location_counter_1;
    end
end

 
%% Make Matrixes for positives and all interactions
%Non-redunent_list_of_all_major_protein_IDs
Protein_MajorProtinID_vert= vertcat(Protein.MajorID_NoIsoforms{:});
Unique_MajorProteinIDs=unique(Protein_MajorProtinID_vert);

%Use all MajorProteinIDs to ensure no protein interactions are miss
Pos_Corum_proteins=Unique_Corum(ismember(Unique_Corum, Unique_MajorProteinIDs)); %Corum part of interactions

%Define variable to write position of proteingroups which are in corum into
count1=1;
length1=length(Corum_Protein_names);
dimension_Protein_MajorID_NoIsoforms=size(Protein.MajorID_NoIsoforms);
Pos=cell(1,1);

%% Find proteins in corum, which are present in dataset
%Create a file to store data
name_processing_file=strcat('Processing data replicate_',mat2str(replicate_counter));
mkdir(name_processing_file);

parfor i=1:length1 %Write out values as i
    Prot1=Corum_Protein_names(i,1);
    Prot2=Corum_Protein_names(i,2);
    AK=(ismember(Prot1, Pos_Corum_proteins)==1)&...
    (ismember(Prot2, Pos_Corum_proteins)==1);
    if AK==1 
        %Map Proteins that are in Corum to Protein.Isoform matrix for Prot1
        for Matrix_check_1_1 =1:dimension_Protein_MajorID_NoIsoforms(1)
        for Matrix_check_1_2 =1:dimension_Protein_MajorID_NoIsoforms(2)
            if isequal(Prot1, Protein.MajorID_NoIsoforms{Matrix_check_1_1,Matrix_check_1_2});
            %Map Proteins that are in Corum to Protein.Isoform matrix for Prot2
            for Matrix_check_2_1 =1:dimension_Protein_MajorID_NoIsoforms(1)
            for Matrix_check_2_2 =1:dimension_Protein_MajorID_NoIsoforms(2)
                if isequal(Prot2, Protein.MajorID_NoIsoforms{Matrix_check_2_1,Matrix_check_2_2});
                    Write_out1=Protein.MajorID_NoIsoforms{Matrix_check_1_1,Matrix_check_1_2};
                    Write_out2=Protein.MajorID_NoIsoforms{Matrix_check_2_1,Matrix_check_2_2};
                %Write out data file
                fileName = strcat(mat2str(i),'_',mat2str(Matrix_check_1_1),'_',mat2str(Matrix_check_2_1),'_Interaction','.csv');
                fid1A = fopen(fileName,'at');
                fprintf(fid1A,'%s,%s,%s,%6.4f,%s,%s,%6.4f,',strcat(Prot1{1},'-',Prot2{1}),...
                    Protein.Isoform{(Matrix_check_1_1(1)),1}, Write_out1{1} ,...
                    Matrix_check_1_1(1),Protein.Isoform{(Matrix_check_2_1(1)),1},...
                    Write_out2{1}, Matrix_check_2_1(1));
                fclose(fid1A);
                movefile(fileName, name_processing_file);
                end
            end
            end
            end
        end
        end
    end
end


%read dictory and concatinate into a single file
%% OutputGaus files filtered out
dirData1 = dir(name_processing_file);      %# Get directory contents
dirData1 = dirData1(~[dirData1.isdir]);  %# Use only the file data
Incorum_fileNames = {dirData1.name};           %# Get file names

ExpName1 = strcat('^(\d+)\_(\d+)\_(\d+)\','_Interaction','.csv');
Incorum_fileData = regexp(Incorum_fileNames,ExpName1,'tokens');  %# Find tokens
index1 = ~cellfun('isempty',Incorum_fileData);        %# Find index of matches
Incorum_fileNames = Incorum_fileNames(index1);        %# remove empty cells
Incorum_fileData = [Incorum_fileData{index1}];                %# Remove non-matches
Incorum_fileData = vertcat(Incorum_fileData{:});             %# Format token data
              
%sorting protein interactions in Corum in correct numerical order
nFiles1 = size(Incorum_fileData,1);              %# Number of files matching format
Incorum_fileData = str2double(Incorum_fileData);        %# Convert from strings to numbers
[Incorum_fileData,index1] = sortrows(Incorum_fileData,1:2);  %# Sort numeric values
Incorum_fileNames =  Incorum_fileNames(index1);               %# Apply sort to file names

Number_of_Protein_interactions_in_corum = length(Incorum_fileNames);
Combined_Incorum_data=cell(Number_of_Protein_interactions_in_corum,7);


for iii = 1:Number_of_Protein_interactions_in_corum 
    %load file of interest
    fileName6 = strcat(name_processing_file,'\',Incorum_fileNames(iii));
    File_to_be_read = fileName6{1};
    
    %import data files
    f = fopen (File_to_be_read);
    Processed_data = textscan(f, '%s %s %s %s %s %s %s', 'Delimiter', ',');
    fclose(f);  
    
    %Copy values to arrays to use in downstream analysis
    %Copy Pos entry
    Pos_entry=Processed_data{1};
    Corum_Dataset{iii}=Pos_entry{1};
    %Copy Protein 1
    prot1_protein_names_entry=Processed_data{2};
    Corum_Dataset_prot1_protein_names{iii}=prot1_protein_names_entry{1};
    %Copy MajorID Protein 1
    prot1_major_names_entry=Processed_data{3};
    Corum_Dataset_prot1_MajorID_names{iii}=prot1_major_names_entry{1};
    %Copy Protein 1 row
    prot1_row_entry=Processed_data{4};
    Corum_Dataset_prot1_row_number(iii)=str2double(prot1_row_entry(1));
    %Copy Protein 2
    prot2_protein_names_entry=Processed_data{5};
    Corum_Dataset_prot2_protein_names{iii}=prot2_protein_names_entry{1};
    %Copy MajorID Protein 2
    prot2_major_names_entry=Processed_data{6};
    Corum_Dataset_prot2_MajorID_names{iii}=prot2_major_names_entry{1};
    %Copy Protein 2 row
    prot2_row_entry=Processed_data{7};
    Corum_Dataset_prot2_row_number(iii)=str2double(prot2_row_entry(1));
    %add Gaus index number
    Known_corum_interaction_index_number(iii)= iii;    
end


%Write out possible interaction list
cd('ROC analysis matrixes');
proteins_interactions_in_corum_name=strcat('proteins_interactions_in_corum replicate_',mat2str(replicate_counter),'.csv');
fid_Eu= fopen(proteins_interactions_in_corum_name,'wt'); % Eu dis interaction output 
fprintf (fid_Eu,'%s,%s,%s,%s,%s,%s,%s,%s,\n','Index interaction number',... %header for output
                'Protein interaction', 'Protein 1 name (with isoform information)', 'Major Protein ID matched in corum of protein 1 (no isoform information)',...
                'Row number of protein 1','Protein 2 name (with isoform information)', 'Major Protein ID matched in corum of protein 2 (no isoform information)',...
                'Row number of protein 2' ); %Write Header
for write_out_loop=1:Number_of_Protein_interactions_in_corum   
fprintf(fid_Eu,'%6f,%s,%s,%s,%6f,%s,%s,%6f,\n',...
    Known_corum_interaction_index_number(write_out_loop),Corum_Dataset{write_out_loop},...
    Corum_Dataset_prot1_protein_names{write_out_loop},...
    Corum_Dataset_prot1_MajorID_names{write_out_loop},...
    Corum_Dataset_prot1_row_number(write_out_loop),...
    Corum_Dataset_prot2_protein_names{write_out_loop},...
    Corum_Dataset_prot2_MajorID_names{write_out_loop},...
    Corum_Dataset_prot2_row_number(write_out_loop));
end            
fclose(fid_Eu);
cd('../');

%% Determine which interactions are unique
%create list with interactions based on row of protein group
Protein_group_interactions.Row_number_grouping=cell(1,1);

for Unique_protein_group_loop=1:Number_of_Protein_interactions_in_corum
    Protein_group_interactions.Row_number_grouping{Unique_protein_group_loop,1}= strcat(num2str(Corum_Dataset_prot1_row_number(1,Unique_protein_group_loop))...
        ,'-',num2str(Corum_Dataset_prot2_row_number(1,Unique_protein_group_loop)));
end

Known_interaction_size=size(Protein_group_interactions.Row_number_grouping);

for  Unique_protein_group_loop=1:Number_of_Protein_interactions_in_corum
    
     %Reset varibles
     Protein_group_uniprot_accessions_interaction=cell(1,1);
     Prot1_Uniprot_accessions=cell(1,1);
     Prot2_Uniprot_accessions=cell(1,1);
     
     %using the major protein groups find the location in corum database array of known interaction within the replicate
     Major_proteingroup_interaction_known_in_corum=Protein_group_interactions.Row_number_grouping{Unique_protein_group_loop,1};
    [Location_of_interaction_in_corum1]=ind2sub(Known_interaction_size, strmatch(Major_proteingroup_interaction_known_in_corum, Protein_group_interactions.Row_number_grouping, 'exact'));
     
     %look up in protein groups and record interactions
     for Unique_protein_group_loop2=1:length(Location_of_interaction_in_corum1);
         Protein_group_interactions.Uniprot_accessions(Unique_protein_group_loop,Unique_protein_group_loop2)=unique(Corum_Dataset(Location_of_interaction_in_corum1(Unique_protein_group_loop2)));
         Protein_group_uniprot_accessions_interaction(Unique_protein_group_loop2)=unique(Corum_Dataset(Location_of_interaction_in_corum1(Unique_protein_group_loop2)));
         Prot1_Uniprot_accessions(1,Unique_protein_group_loop2)=unique(Corum_Dataset_prot1_MajorID_names(Location_of_interaction_in_corum1(Unique_protein_group_loop2)));
         Prot2_Uniprot_accessions(1,Unique_protein_group_loop2)=unique(Corum_Dataset_prot2_MajorID_names(Location_of_interaction_in_corum1(Unique_protein_group_loop2))); 
     end
     
     %Unique uniprot accession numbers
     unique_Uniprot_accessions_interaction=unique(Protein_group_uniprot_accessions_interaction);
     lenght_unique_Uniprot_accessions_interaction=length(unique_Uniprot_accessions_interaction);
     for Unique_protein_group_loop3=1:lenght_unique_Uniprot_accessions_interaction
         Protein_group_interactions.Unique_Uniprot(Unique_protein_group_loop,Unique_protein_group_loop3)=unique_Uniprot_accessions_interaction(Unique_protein_group_loop3);
     end
     
     %record unique Prot 1 uniprot entries
     unique_Prot1_Uniprot_accessions=unique(Prot1_Uniprot_accessions);
     lenght_unique_Prot1_Uniprot_accessions=length(unique_Prot1_Uniprot_accessions);
     for Unique_protein_group_loop4=1:lenght_unique_Prot1_Uniprot_accessions
         Protein_group_interactions.Uniprot_Prot1(Unique_protein_group_loop,Unique_protein_group_loop4)=unique_Prot1_Uniprot_accessions(Unique_protein_group_loop4);
     end
     
     %record unique Prot 2 uniprot entries
     unique_Prot2_Uniprot_accessions=unique(Prot2_Uniprot_accessions);
     lenght_unique_Prot2_Uniprot_accessions=length(unique_Prot2_Uniprot_accessions);
     for Unique_protein_group_loop5=1:lenght_unique_Prot2_Uniprot_accessions
         Protein_group_interactions.Uniprot_Prot2(Unique_protein_group_loop,Unique_protein_group_loop5)=unique_Prot2_Uniprot_accessions(Unique_protein_group_loop5);
     end
     
     %Record interaction entries observed
     Protein_group_interactions.multiple_interactions(Unique_protein_group_loop,lenght_unique_Uniprot_accessions_interaction)=1;
     
end

%Generate headers to write out
size_multiple_interaction=size(Protein_group_interactions.multiple_interactions);
size_uniprot_Prot1=size(Protein_group_interactions.Uniprot_Prot1);
size_uniprot_Prot2=size(Protein_group_interactions.Uniprot_Prot2);
size_unique_uniprot=size(Protein_group_interactions.Unique_Uniprot);

for column_counter=1:size_multiple_interaction(2)
    Column_name1{column_counter}=strcat('Interaction matches majorprotein ID combination_',mat2str(column_counter));
end

for column_counter=1:size_uniprot_Prot1(2)
    Column_name2{column_counter}=strcat('Unique Uniprot accession protein1_',mat2str(column_counter));
end

for column_counter=1:size_uniprot_Prot2(2)
    Column_name3{column_counter}=strcat('Unique Uniprot accession protein2_',mat2str(column_counter));
end

for column_counter=1:size_unique_uniprot(2)
    Column_name4{column_counter}=strcat('Unique Uniprot accession_',mat2str(column_counter));
end


Column_header1=repmat('%s,', 1, (size_unique_uniprot(2)+...
size_uniprot_Prot2(2)+size_uniprot_Prot1(2)+size_multiple_interaction(2)+1));

%Write out possible interaction list
cd('ROC analysis matrixes');
Unique_proteins_interactions_in_corum_name=strcat('Unique_proteins_interactions_in_corum_',mat2str(replicate_counter),'.csv');
fid_Eu= fopen(Unique_proteins_interactions_in_corum_name,'wt'); % Eu dis interaction output 
fprintf (fid_Eu,[Column_header1, '\n'],... %header for output
'Protein Group row number',Column_name4{:},...
 Column_name2{:},Column_name3{:},Column_name1{:}); %Write Header
for write_out_loop2=1:Number_of_Protein_interactions_in_corum   
fprintf(fid_Eu,'%s,',Protein_group_interactions.Row_number_grouping{write_out_loop2,1});
fprintf(fid_Eu,'%s,', Protein_group_interactions.Unique_Uniprot{write_out_loop2,1:end});
fprintf(fid_Eu,'%s,', Protein_group_interactions.Uniprot_Prot1{write_out_loop2,1:end});
fprintf(fid_Eu,'%s,', Protein_group_interactions.Uniprot_Prot2{write_out_loop2,1:end});
fprintf(fid_Eu,'%6.4f,', Protein_group_interactions.multiple_interactions(write_out_loop2,1:end));
fprintf(fid_Eu,',\n');
end            
fclose(fid_Eu);
cd('../');

%Create list of all corum interactions 
%combine values together
dimension_Unique_Uniprot=size(Protein_group_interactions.Unique_Uniprot);
Combined_Protein_group_interactions=zeros(0,0);
for combined_counter=1:dimension_Unique_Uniprot(2)
Combined_Protein_group_interactions=[Combined_Protein_group_interactions;Protein_group_interactions.Unique_Uniprot(:,combined_counter)];
end
Size_of_combined_protein_group_interactions=length(Combined_Protein_group_interactions);
%Create array
Cells_Combined_Protein_group_interactions=cell(1,1);
%Create counter
Combined_Protein_group_interactions_Counter=1;
%Test if cells are empty, if not add to array
for copy_cell_loop_counter=1:Size_of_combined_protein_group_interactions
    if ~isempty(Combined_Protein_group_interactions{copy_cell_loop_counter})
        Cells_Combined_Protein_group_interactions{Combined_Protein_group_interactions_Counter,1}=Combined_Protein_group_interactions{copy_cell_loop_counter};
        Combined_Protein_group_interactions_Counter=Combined_Protein_group_interactions_Counter+1;
    end
end

%Copy unique entries from Protein_group_interactions 
Unique_combined_Protein_group_interactions= unique(Cells_Combined_Protein_group_interactions);

%array with alterative MajorID to check
Number_MajorIDs=zeros(Dimensions_Gaus_import,1);

%Create counter of alterative MajorID to check
for counter_1=1:Dimensions_Gaus_import
    for counter_2=1:Dimension_of_Protein_IDs(2)
Number_MajorIDs(counter_1)=Number_MajorIDs(counter_1)+ ~isempty(Protein.MajorID_NoIsoforms{counter_1,counter_2});
    end
end

% Make positive matrix
TP_Matrix= zeros(Dimensions_Gaus_import,Dimensions_Gaus_import);
% Make matrix of protein group with multiple possible positive interactions
Multiple_known_interaction_in_proteingroup_Matrix= zeros(Dimensions_Gaus_import,Dimensions_Gaus_import);

h = waitbar(0,'Calculating True positive matrix');

for i=1:Dimensions_Gaus_import
    %check if the i counter is n/100 of the Dimensions_Gaus_import, if so change counter
    if (round(i/(Dimensions_Gaus_import/100))==i/(Dimensions_Gaus_import/100))
    waitbar(i/Dimensions_Gaus_import,h,sprintf('%d%% of True positive matrix caliculated',round((i/Dimensions_Gaus_import)*100)));
    end
    %begin actually loop
    for i_counter=1:Number_MajorIDs(i)
        %check if protein is in corum
        if  ismember(Protein.MajorID_NoIsoforms{i,i_counter},Unique_Corum);
        parfor ii=1:Dimensions_Gaus_import   
        for ii_counter=1:Number_MajorIDs(ii)
            if  ismember(Protein.MajorID_NoIsoforms{ii,ii_counter},Unique_Corum);
            %using string is slow, use the row number instead
             Interaction= strcat(Protein.MajorID_NoIsoforms{i,i_counter},'-',Protein.MajorID_NoIsoforms{ii,ii_counter});    
             AK=(ismember(Interaction, Unique_combined_Protein_group_interactions)==1);
                if AK==1 && TP_Matrix(i,ii)==0
                TP_Matrix(i,ii)=1;
                elseif AK==1 && TP_Matrix(i,ii)==1    
                Multiple_known_interaction_in_proteingroup_Matrix(i,ii)=1;
                end
            end
        end
        end
        end
    end
end
close(h);

%write True positive matrix to file to enable file to be reloaded quickly
%for re-analysis

cd('ROC analysis matrixes');
TP_name=strcat('TP_matrix_replicate',mat2str(replicate_counter),'.csv');
Multiple_TP_known_interaction_name=strcat('Multiple_known_TP_interaction_in_proteingroup_Matrix_replicate',mat2str(replicate_counter),'.csv');
csvwrite(TP_name,TP_Matrix);
csvwrite(Multiple_TP_known_interaction_name,Multiple_known_interaction_in_proteingroup_Matrix);
cd('../');

%to reload positive matrix run import code below
%cd('ROC analysis matrixes');
%TP_name=strcat('TP_matrix_replicate',mat2str(replicate_counter),'.csv');
%Multiple_TP_known_interaction_name=strcat('Multiple_known_TP_interaction_in_proteingroup_Matrix_replicate',mat2str(replicate_counter),'.csv');
%TP_Matrix = importdata(TP_name);
%Multiple_known_interaction_in_proteingroup_Matrix =importdata(Multiple_TP_known_interaction_name);
%cd('../');

% Calculate possible interaction matrix, asks are both proteins in corum
Int_matrix= zeros(Dimensions_Gaus_import,Dimensions_Gaus_import);
% Make matrix of protein group with multiple possible positive interactions
Multiple_matches_within_protein_group_Matrix = zeros(Dimensions_Gaus_import,Dimensions_Gaus_import);

for i=1:Dimensions_Gaus_import
    %i %Display protein number the loop is up to to enable users to monitor progress
    for i_counter=1:Number_MajorIDs(i)
    if ismember(Protein.MajorID_NoIsoforms{i,i_counter}, Pos_Corum_proteins)
        parfor ii=1:Dimensions_Gaus_import
        for ii_counter=1:Number_MajorIDs(ii)
            if ismember(Protein.MajorID_NoIsoforms{ii,ii_counter}, Pos_Corum_proteins)
            AK=((ismember(Protein.MajorID_NoIsoforms{i,i_counter}, Pos_Corum_proteins)==1)&...
            (ismember(Protein.MajorID_NoIsoforms{ii,ii_counter}, Pos_Corum_proteins)==1));
                if AK==1 && Int_matrix(i,ii)== 0
                Int_matrix(i,ii)=1;
                elseif AK==1 && Int_matrix(i,ii)== 1 
                Multiple_matches_within_protein_group_Matrix(i,ii)=1;
                end
            end
        end
        end
    end
    end
end

%write possible interaction matrix to file to enable file to be quickly reloaded
%for re-analysis
    
cd('ROC analysis matrixes');
INT_matrix_name=strcat('INT_matrix_replicate',mat2str(replicate_counter),'.csv');
Multiple_known_INT_interaction_name=strcat('Multiple_known_INT_interaction_in_proteingroup_Matrix_replicate',mat2str(replicate_counter),'.csv');
csvwrite(INT_matrix_name,Int_matrix);
csvwrite(Multiple_known_INT_interaction_name,Multiple_matches_within_protein_group_Matrix);
cd('../');


%to reload positive matrix run import code below
%cd('ROC analysis matrixes');
%INT_matrix_name=strcat('INT_matrix_replicate',mat2str(replicate_counter),'.csv');
%Multiple_known_INT_interaction_name=strcat('Multiple_known_INT_interaction_in_proteingroup_Matrix_replicate',mat2str(replicate_counter),'.csv');
%Int_matrix = importdata(INT_matrix_name);
%Multiple_matches_within_protein_group_Matrix=importdata(Multiple_known_INT_interaction_name);
%cd('../');

% Create self interaction matrix
Self_Int_matrix= zeros(Dimensions_Gaus_import,Dimensions_Gaus_import);
for i=1:Dimensions_Gaus_import
    parfor ii=1:Dimensions_Gaus_import
        Protein.NoIsoform;
        NS=isequal(Protein.NoIsoform{i},Protein.NoIsoform{ii});
        if NS==1
            Self_Int_matrix(i,ii)=1;
        end

    end
end
%write possible interaction matrix to file to enable file to be quickly reloaded
%for re-analysis

cd('ROC analysis matrixes');
Self_interaction_matrix_name=strcat('Self_interaction_matrix_replicate',mat2str(replicate_counter),'.csv');
csvwrite(Self_interaction_matrix_name,Self_Int_matrix);
cd('../');

%to reload positive matrix run import code below
%cd('ROC analysis matrixes');
%Self_interaction_matrix_name=strcat('Self_interaction_matrix_replicate',mat2str(replicate_counter),'.csv');
%Self_Int_matrix = importdata(Self_interaction_matrix_name);
%cd('../');

% Make false positive matrix, proteins that are both in corum but no
% interaction known
Inverse_TP_matrix=~TP_Matrix;
FP_Matrix= (Int_matrix & Inverse_TP_matrix);

%write False positive matrix to file to enable file to be reloaded quickly
%for re-analysis

cd('ROC analysis matrixes');
FP_matrix_name=strcat('FP_matrix_replicate',mat2str(replicate_counter),'.csv');
csvwrite(FP_matrix_name,FP_Matrix);
cd('../');

%to reload positive matrix run import code below
%cd('ROC analysis matrixes');
%FP_Matrix = importdata('FP_matrix.csv');
%cd('../');

%Possible interaction
inverse_self=~Self_Int_matrix; %creation of inverse of self interactions

%% Write out statistic
cd('ROC analysis matrixes');
%create list of unique MajorIDs with isoform information
Protein_MajorProtinIDwith_isoforms_vert= vertcat(Protein.MajorIDs{:});
Unique_MajorProteinIDswith_isoforms=unique(Protein_MajorProtinIDwith_isoforms_vert);

%count the number of proteins/interactions in corum
number_of_Unique_Corum=length(Unique_Corum); % number of unique entires in corum
number_of_Unique_MajorProteinIDs=length(Unique_MajorProteinIDs); % number of unique MajorProteinIDs (No isoform information) within dataset
number_of_Unique_MajorProteinIDswith_isoforms=length(Unique_MajorProteinIDswith_isoforms); % number of unique MajorProteinIDs (with isoform information) within dataset

number_of_pos_corum_proteins=length(Pos_Corum_proteins); % number of entries in corum also within dataset
number_of_Cells_Combined_Protein_group_interactions=length(Cells_Combined_Protein_group_interactions); %redundant list of all interactions in corum possible with protei Identifications
number_of_Unique_combined_Protein_group_interactions=length(Unique_combined_Protein_group_interactions); %non-redundant list of all interactions in corum possible with protei Identifications

%count interactions from matrix
inverse_self_triu=triu(inverse_self);
number_of_possible_interactions= length(find(inverse_self_triu==1));
output_matrix_01 = (Int_matrix & inverse_self);
output_matrix_01_triu = triu(output_matrix_01);
number_of_possible_interactions_of_corum_proteins= length(find(output_matrix_01_triu==1));
output_matrix_02 = (Multiple_known_interaction_in_proteingroup_Matrix & inverse_self);
output_matrix_02_triu = triu(output_matrix_02);
number_of_Multiple_known_interaction_in_proteingroup_Matrix= length(find(output_matrix_02_triu==1));
output_matrix_03 = (Multiple_matches_within_protein_group_Matrix & inverse_self);
output_matrix_03_triu = triu(output_matrix_03);
number_of_Multiple_matches_within_protein_group_Matrix= length(find(output_matrix_03_triu==1));
number_of_self_interactions= length(find(Self_Int_matrix==1));
output_matrix_04 = (TP_Matrix & inverse_self);
output_matrix_04_triu = triu(output_matrix_04);
number_of_TP_interactions= length(find(output_matrix_04_triu==1));
output_matrix_05 = (FP_Matrix & inverse_self);
output_matrix_05_triu = triu(output_matrix_05);
number_of_TN_interactions= length(find(output_matrix_05_triu==1));

%Write to file
Corum_vs_dataset_statistic_name=strcat('Corum_vs_dataset_statistic_replicate',mat2str(replicate_counter),'.csv');
fid_stats= fopen(Corum_vs_dataset_statistic_name,'wt'); 
fprintf (fid_stats,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',... %header for output
                'Number of unique entries in corum',....
                'Number of unique MajorProteinIDs with isoform information in dataset',...
                'Number of unique MajorProteinIDs without isoform information in dataset',...
                'Unique Corum entries within dataset', 'Redundant Corum interactions within dataset',...
                'Unique Corum interactions within dataset','Number of possible interactions',...
                'Number of possible interactions of proteins in corum','Maximum number of TP',...
                'Maximum number of TN', 'Number of self interactions',...
                'Number of proteingroups with members that have mutliple known interactions in corum',...
                'Number of proteingroups with members that have multiple interactions in corum'); %Write Header           
fprintf(fid_stats,'%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',...
    number_of_Unique_Corum(1), number_of_Unique_MajorProteinIDswith_isoforms(1),...
    number_of_Unique_MajorProteinIDs(1), number_of_pos_corum_proteins(1),...
    number_of_Cells_Combined_Protein_group_interactions(1),number_of_Unique_combined_Protein_group_interactions(1),...
    number_of_possible_interactions(1), number_of_possible_interactions_of_corum_proteins(1) , number_of_TP_interactions(1),...
    number_of_TN_interactions(1), number_of_self_interactions(1),number_of_Multiple_known_interaction_in_proteingroup_Matrix(1),... 
    number_of_Multiple_matches_within_protein_group_Matrix(1));           
fclose(fid_stats);
cd('../');

% Save figure
cd('ROC analysis related Figures');

%% PCA analysis (This module is designed to provide users an alterative means to assass there data) 

%Analysis each variable with PCA
PCA_analysis_Eu =myPCA(Chromatograms', 2);
PCA_analysis_normalsied_gaus =myPCA(gaussian_fit', 2);

%Write out PCA analysis
PCA_txt_name1=strcat('PCA_Chromatograms_analysis_replicates_',mat2str(replicate_counter),'.csv');
PCA_list1= fopen(PCA_txt_name1,'wt'); % create the output file with the header infromation
fprintf (PCA_list1,'%s,%s,%s,%s,%s\n',... %header for OutputGaus output
                'Protein name', 'Unique_identifier','Protein Center','PCA component 1',...
                'PCA component 2'); %Write Header                     
for writeout_counter1= 1:Len(1)
fprintf (PCA_list1,'%s,%s,%6.4f,%6.4f,%6.4f\n',...   
Protein.Isoform{writeout_counter1,1},Protein.Unique_ID{writeout_counter1,1},...
Protein.Center(writeout_counter1,1),PCA_analysis_Eu(1,writeout_counter1),...
PCA_analysis_Eu(2,writeout_counter1));
end
fclose(PCA_list1);

%Caliculate dimensions for density analysis
PCA_dimensionA1=round(abs(min(PCA_analysis_Eu(1,:)))+ abs(max(PCA_analysis_Eu(1,:))));
PCA_dimensionA2=round(abs(min(PCA_analysis_Eu(2,:)))+ abs(max(PCA_analysis_Eu(2,:))));
%generate density measurement of PCA anlaysis

xi=linspace(min(PCA_analysis_Eu(1,:)),max(PCA_analysis_Eu(1,:)),100);
yi=linspace(min(PCA_analysis_Eu(2,:)),max(PCA_analysis_Eu(2,:)),100);

%yi corresponds to the replicate being adjusted
%xi corresponds to the replicate being used as the reference
xr=interp1(xi, 1:numel(xi),PCA_analysis_Eu(1,:),'nearest');
yr=interp1(yi, 1:numel(yi),PCA_analysis_Eu(2,:),'nearest');

Z=accumarray([xr' yr'], 1, [100 100]);

%Convert zeros to NaN
Z(Z==0)=NaN;

figure
f_pca1=subplot(2,1,1)
scatter(PCA_analysis_Eu(2,:), PCA_analysis_Eu(1,:),10,...
              'MarkerFaceColor','b',...
              'LineWidth',1);
title('PCA analysis -Chromatograms (Scatter plot)','FontSize', 12);
ylabel('PCA component 1','FontSize', 8);
xlabel('PCA component 2','FontSize', 8);

f_pca1=subplot(2,1,2)
contourf(Z);
colormap(summer);
caxis([1 10]);
title('PCA analysis -Chromatograms (Density map)','FontSize', 12);
ylabel('PCA component 1','FontSize', 8);
xlabel('PCA component 2','FontSize', 8);
colormap summer
colorbar('YTickLabel',...
    {'2','4','6','8','>=10'},...
     'location','eastoutside');

PCA_figure=strcat('PCA_Chromatograms_analysis_replicates_',mat2str(replicate_counter),'.png');
saveas(f_pca1, PCA_figure);


%Write out PCA analysis
PCA_txt_name2=strcat('PCA_Normalized_Gaus_analysis_replicates_',mat2str(replicate_counter),'.csv');
PCA_list2= fopen(PCA_txt_name2,'wt'); % create the output file with the header infromation
fprintf (PCA_list2,'%s,%s,%s,%s,%s\n',... %header for OutputGaus output
                'Protein name', 'Unique_identifier','Protein Center','PCA component 1',...
                'PCA component 2'); %Write Header                     
for writeout_counter1= 1:Len(1)
fprintf (PCA_list2,'%s,%s,%6.4f,%6.4f,%6.4f\n',...   
Protein.Isoform{writeout_counter1,1},Protein.Unique_ID{writeout_counter1,1},...
Protein.Center(writeout_counter1,1),...
PCA_analysis_normalsied_gaus(1,writeout_counter1),...
PCA_analysis_normalsied_gaus(2,writeout_counter1));
end
fclose(PCA_list2);

%Caliculate dimensions for density analysis
PCA_dimensionA1=round(abs(min(PCA_analysis_normalsied_gaus(1,:)))+ abs(max(PCA_analysis_normalsied_gaus(1,:))));
PCA_dimensionA2=round(abs(min(PCA_analysis_normalsied_gaus(2,:)))+ abs(max(PCA_analysis_normalsied_gaus(2,:))));
%generate density measurement of PCA anlaysis

xi=linspace(min(PCA_analysis_normalsied_gaus(1,:)),max(PCA_analysis_normalsied_gaus(1,:)),100);
yi=linspace(min(PCA_analysis_normalsied_gaus(2,:)),max(PCA_analysis_normalsied_gaus(2,:)),100);

%yi corresponds to the replicate being adjusted
%xi corresponds to the replicate being used as the reference
xr=interp1(xi, 1:numel(xi),PCA_analysis_normalsied_gaus(1,:),'nearest');
yr=interp1(yi, 1:numel(yi),PCA_analysis_normalsied_gaus(2,:),'nearest');

Z=accumarray([xr' yr'], 1, [100 100]);

%Convert zeros to NaN
Z(Z==0)=NaN;

figure
f_pca2=subplot(2,1,1)
scatter(PCA_analysis_normalsied_gaus(2,:), PCA_analysis_normalsied_gaus(1,:),10,...
              'MarkerFaceColor','b',...
              'LineWidth',1);
title('PCA analysis -normalized Gaussian profile (Scatter plot)','FontSize', 12);
ylabel('PCA component 1','FontSize', 8);
xlabel('PCA component 2','FontSize', 8);

f_pca2=subplot(2,1,2)
contourf(Z);
colormap(summer);
caxis([1 10]);
title('PCA analysis -normalized Gaussian profile (Density map)','FontSize', 12);
ylabel('PCA component 1','FontSize', 8);
xlabel('PCA component 2','FontSize', 8);
colormap summer
colorbar('YTickLabel',...
    {'2','4','6','8','>=10'},...
     'location','eastoutside');

PCA_figure=strcat('PCA_Normalized_Gaussian_analysis_replicates_',mat2str(replicate_counter),'.png');
saveas(f_pca2, PCA_figure);

% format array into correct format
PCA_analysis_Eu=PCA_analysis_Eu';
PCA_analysis_normalsied_gaus=PCA_analysis_normalsied_gaus';

% calculate similarity
Dist.PCA_Chromatograms = squareform(pdist(PCA_analysis_Eu,'euclidean'));
Dist.PCA_normalized_gaussian = squareform(pdist(PCA_analysis_normalsied_gaus,'euclidean'));

%% TSNE analysis (This module is designed to provide users an alterative means to assass there data) 

%Analysis each variable with TSNE
tsne_analysis_Eu =tsne(Chromatograms);
tsne_analysis_normalsied_gaus =tsne(gaussian_fit);

%Write out T-SNE analysis
tsne_txt_name1=strcat('T-SNE_Chromatograms_analysis_replicates_',mat2str(replicate_counter),'.csv');
tsne_list1= fopen(tsne_txt_name1,'wt'); % create the output file with the header infromation
fprintf (tsne_list1,'%s,%s,%s,%s,%s\n',... %header for OutputGaus output
                'Protein name', 'Unique_identifier','Protein Center','T-SNE component 1',...
                'T-SNE component 2'); %Write Header                     
for writeout_counter1= 1:Len(1)
fprintf (tsne_list1,'%s,%s,%6.4f,%6.4f,%6.4f\n',...   
Protein.Isoform{writeout_counter1,1},Protein.Unique_ID{writeout_counter1,1},...
Protein.Center(writeout_counter1,1),tsne_analysis_Eu(writeout_counter1,1),...
tsne_analysis_Eu(writeout_counter1,2));
end
fclose(tsne_list1);

%Caliculate dimensions for density analysis
tsne_dimensionA1=round(abs(min(tsne_analysis_Eu(:,1)))+ abs(max(tsne_analysis_Eu(:,1))));
tsne_dimensionA2=round(abs(min(tsne_analysis_Eu(:,2)))+ abs(max(tsne_analysis_Eu(:,2))));

%generate density measurement of PCA anlaysis
xi=linspace(min(tsne_analysis_Eu(:,1)),max(tsne_analysis_Eu(:,1)),tsne_dimensionA1);
yi=linspace(min(tsne_analysis_Eu(:,2)),max(tsne_analysis_Eu(:,2)),tsne_dimensionA2);

%yi corresponds to the replicate being adjusted
%xi corresponds to the replicate being used as the reference
xr=interp1(xi, 1:numel(xi),tsne_analysis_Eu(:,1),'nearest');
yr=interp1(yi, 1:numel(yi),tsne_analysis_Eu(:,2),'nearest');

Z=accumarray([xr yr], 1, [tsne_dimensionA1 tsne_dimensionA2]);

%Convert zeros to NaN
Z(Z==0)=NaN;

figure
f_tsne1=subplot(2,1,1)
scatter(tsne_analysis_Eu(:,2), tsne_analysis_Eu(:,1),10,...
              'MarkerFaceColor','b',...
              'LineWidth',1);
title('T-SNE analysis -Chromatograms (Scatter plot)','FontSize', 12);
ylabel('T-SNE component 1','FontSize', 8);
xlabel('T-SNE component 2','FontSize', 8);

f_tsne1=subplot(2,1,2)
contourf(Z);
colormap(summer);
caxis([1 10]);
title('T-SNE analysis -Chromatograms (Density map)','FontSize', 12);
ylabel('T-SNE component 1','FontSize', 8);
xlabel('T-SNE component 2','FontSize', 8);
ylim([-10 tsne_dimensionA1*1.1]);
xlim([-10 tsne_dimensionA2*1.1]);
colormap summer
colorbar('YTickLabel',...
    {'2','4','6','8','>=10'},...
     'location','eastoutside');

TSNE_figure=strcat('TSNE_Chromatograms_analysis_replicates_',mat2str(replicate_counter),'.png');
saveas(f_tsne1, TSNE_figure);


%Write out T-SNE analysis
tsne_txt_name2=strcat('T-SNE_Normalsied_gaus_analysis_replicates_',mat2str(replicate_counter),'.csv');
tsne_list2= fopen(tsne_txt_name2,'wt'); % create the output file with the header infromation
fprintf (tsne_list2,'%s,%s,%s,%s,%s\n',... %header for OutputGaus output
                'Protein name', 'Unique_identifier','Protein Center','T-SNE component 1',...
                'T-SNE component 2'); %Write Header                     
for writeout_counter1= 1:Len(1)
fprintf (tsne_list2,'%s,%s,%6.4f,%6.4f,%6.4f\n',...   
Protein.Isoform{writeout_counter1,1},Protein.Unique_ID{writeout_counter1,1},...
Protein.Center(writeout_counter1,1),tsne_analysis_normalsied_gaus(writeout_counter1,1),...
tsne_analysis_normalsied_gaus(writeout_counter1,2));
end
fclose(tsne_list2);

%Caliculate dimensions for density analysis
tsne_dimensionA1=round(abs(min(tsne_analysis_normalsied_gaus(:,1)))+ abs(max(tsne_analysis_normalsied_gaus(:,1))));
tsne_dimensionA2=round(abs(min(tsne_analysis_normalsied_gaus(:,2)))+ abs(max(tsne_analysis_normalsied_gaus(:,2))));

%generate density measurement of PCA anlaysis
xi=linspace(min(tsne_analysis_normalsied_gaus(:,1)),max(tsne_analysis_normalsied_gaus(:,1)),tsne_dimensionA1);
yi=linspace(min(tsne_analysis_normalsied_gaus(:,2)),max(tsne_analysis_normalsied_gaus(:,2)),tsne_dimensionA2);

%yi corresponds to the replicate being adjusted
%xi corresponds to the replicate being used as the reference
xr=interp1(xi, 1:numel(xi),tsne_analysis_normalsied_gaus(:,1),'nearest');
yr=interp1(yi, 1:numel(yi),tsne_analysis_normalsied_gaus(:,2),'nearest');

Z=accumarray([xr yr], 1, [tsne_dimensionA1 tsne_dimensionA2]);

%Convert zeros to NaN
Z(Z==0)=NaN;

figure
f_tsne2=subplot(2,1,1)
scatter(tsne_analysis_normalsied_gaus(:,2), tsne_analysis_normalsied_gaus(:,1),10,...
              'MarkerFaceColor','b',...
              'LineWidth',1);
title('T-SNE analysis -Normalsied Gaussian (Scatter plot)','FontSize', 12);
ylabel('T-SNE component 1','FontSize', 8);
xlabel('T-SNE component 2','FontSize', 8);

f_tsne2=subplot(2,1,2)
contourf(Z);
colormap(summer);
caxis([1 10]);
title('T-SNE analysis -Normalsied Gaussian (Density map)','FontSize', 12);
ylabel('T-SNE component 1','FontSize', 8);
xlabel('T-SNE component 2','FontSize', 8);
ylim([-10 tsne_dimensionA1*1.1]);
xlim([-10 tsne_dimensionA2*1.1]);
colormap summer
colorbar('YTickLabel',...
    {'2','4','6','8','>=10'},...
     'location','eastoutside');

TSNE_figure=strcat('TSNE_Normalsied_Gaussian_analysis_replicates_',mat2str(replicate_counter),'.png');
saveas(f_tsne2, TSNE_figure);


% calculate similarity
Dist.tsne_Chromatograms = squareform(pdist(tsne_analysis_Eu,'euclidean'));
Dist.tsne_normalized_gaussian = squareform(pdist(tsne_analysis_normalsied_gaus,'euclidean'));


%% ROC analysis
%Loop to determine FPR TPR precision recall
fn = fieldnames(Dist);

for i=1:length(fn)
    
count=0;
nr=500;
Recall_out.(fn{i})=zeros(nr,1);
Precision_out.(fn{i})=zeros(nr,1);
TPR_out.(fn{i})=zeros(nr,1);
FPR_out.(fn{i})=zeros(nr,1);

filename = strcat((fn{i}),'_replicate',mat2str(replicate_counter),'.csv');
fid5 = fopen(filename,'w');
fprintf (fid5,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n',...
        'Distance', 'TP', 'FP', 'FN', 'TN', 'Interactions within defined distance','Interactions within defined distance also in Corum','Recall', 'Precision', 'FPR', 'TPR'); %Write header


    for Distance= max(max(Dist.(fn{i})))/1000:max(max(Dist.(fn{i})))/1000:max(max(Dist.(fn{i}))) %Loop to calculate performance of variables  %for Distance= 0.1:0.1:50 %Loop to calculate performance of variables       
        MatrixGaussian = Dist.(fn{i}) < Distance;
        count=1+count;
        
        %As Euclidean distance represents a measurement across the total
        %graident filter out measurements of guassian< 2 fractions away
        if isequal(fn(i), fn(1))  
        MatrixGaussian= (MatrixGaussian & (Dist.(fn{2}) < 2));
        
        else
        MatrixGaussian = Dist.(fn{i}) < Distance;
        
        end
        
        %Total number of interaction within distance
        Total_Int = (MatrixGaussian & inverse_self);
        Total_Int2_Triu =triu(Total_Int); %Only consider interactions in the top right of the matrix
        Total_Int_number= length(find(Total_Int2_Triu==1));
        
        %Protein interactions of protein in Corum within distance
        Possible_Int1 = (MatrixGaussian & Int_matrix);
        Possible_Int2 = (Possible_Int1 & inverse_self);
        Possible_Int2_Triu = triu(Possible_Int2); %Only consider interactions in the top right of the matrix
        Possible_Int_number= length(find(Possible_Int2_Triu==1));
        
        % Calculate True possitives
        TP_Matrix_Int1 = (Possible_Int1 & TP_Matrix);
        TP_Matrix_Int2 = (TP_Matrix_Int1 & inverse_self);
        TP_Matrix_Int2_Triu = triu(TP_Matrix_Int2); %Only consider interactions in the top right of the matrix
        TP= length(find(TP_Matrix_Int2_Triu==1));

        % Calculate false negatives
        Inverse_TP_matrix1=~TP_Matrix_Int2;
        MaxTP_matrix1 = (TP_Matrix & Inverse_TP_matrix1);
        MaxTP_matrix2 = (MaxTP_matrix1 & inverse_self);    
        MaxTP_matrix2_Triu = triu(MaxTP_matrix2); %Only consider interactions in the top right of the matrix
        FN= length(find(MaxTP_matrix2_Triu==1));
        
        % Calculate false positives 
        FP_Matrix_Int1= (Possible_Int2 & FP_Matrix);
        FP_Matrix_Int2= (FP_Matrix_Int1 & inverse_self);
        FP_Matrix_Int2_Triu = triu(FP_Matrix_Int2); %Only consider interactions in the top right of the matrix        
        FP= length(find(FP_Matrix_Int2_Triu==1));
        
        % Find True negatives
        Inverse_false_postive_matrix=~FP_Matrix_Int1;
        Int_matrix_minus_TP = (Int_matrix & Inverse_TP_matrix);
        Int5_matrix= (Int_matrix_minus_TP & Inverse_false_postive_matrix);
        Int6_matrix = (Int5_matrix & inverse_self);
        Int6_matrix_Triu = triu(Int6_matrix); %Only consider interactions in the top right of the matrix 
        TN= (length(find(Int6_matrix_Triu==1)));

        Recall = TP/(TP+FN);
        Precision = TP/(TP+FP);
        TPR = TP/(TP+FN);
        FPR =FP/(FP+TN);

        Recall_out.(fn{i})(count,1)= Recall;
        Precision_out.(fn{i})(count,1)=Precision;
        TPR_out.(fn{i})(count,1)=TPR;
        FPR_out.(fn{i})(count,1)=FPR;
        
        %Write to file        
        
        fprintf (fid5,'%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f\n',...
            Distance, TP, FP, FN, TN, Total_Int_number,Possible_Int_number, Recall, Precision, FPR, TPR);
        
    end
    fclose(fid5);
    
    %% Make plots
    filename1 = strcat((fn{i}),'_replicate',mat2str(replicate_counter));
    Figure.(fn{i}) = figure('Color',[1 1 1],'Visible','off');    
    subplot(2,2,1)
    plot(Recall_out.(fn{i}),Precision_out.(fn{i}),'linewidth',2, 'color',[0,0,0]);
    title(fn{i})
    ylabel({'Precision'},'FontSize',10);
    xlabel({'Recall'},'FontSize',10);
    set(gca,'XTick',[0.00 0.25 0.50 0.75 1.00])
    

    subplot(2,2,2)
    plot(FPR_out.(fn{i}),TPR_out.(fn{i}),'linewidth',2, 'color',[0,0,0]);
    line([0,1],[0,1],'linewidth',1, 'color',[0,0,0])
    title(fn{i})
    ylabel({'True positive rate (sensitivity)'},'FontSize',10);
    xlabel({'False positive rate (1-specificity)'},'FontSize',10);
    axis([0,1,0,1])
    set(gca,'YTick',[0.00 0.25 0.50 0.75 1.00])
    set(gca,'XTick',[0.00 0.25 0.50 0.75 1.00])
    print(Figure.(fn{i}),'-dpdf','-r600', filename1)

end
cd('../');

%Fitting Eucledian distance curve to find the value which gives a precision
%of 0.8
import_name=strcat('Euc_replicate',mat2str(replicate_counter),'.csv');

%Import results from initial analysis
cd('ROC analysis related Figures'); % move into the ROC figure folder to import Euc file
Eucldian_distance_import = importdata(import_name);%import Eucldian distance values
cd('../');

cd('ROC analysis data');

%Test if the Euclidian distance ever gives a precision above 0.8
if max(Eucldian_distance_import.data(:,9)>0.80)
[~, idx] = min(abs(Eucldian_distance_import.data(:,9) - 0.8)); %determines the postion which gives a precision closest to 0.80
else
idx=1; % set the distance to the tightest setting as the 
end

%Determine interactions within precision of >0.80 based on Eu distance
Eucldian_distance_with_precision_abover_0_8= Eucldian_distance_import.data(idx,1); %determines the Eucldian distance which gives a precision closest to 0.8
Interactions_detected_due_to_High_Eucldian_distance = Dist.Euc <= Eucldian_distance_with_precision_abover_0_8;    
Interactions_Based_only_Eucldian_distance = ((Interactions_detected_due_to_High_Eucldian_distance & inverse_self) & (Dist.Center <= 2));

%create matrix for Eu dis interactions with precision >0.8
Final_Interactions_Based_only_Eucldian_distance=zeros(Dimensions_Gaus_import,Dimensions_Gaus_import);
Eu_distance_matched=cell(1,1);

internal_counter=1;
for Eucldian_counter1=1: Dimensions_Gaus_import % need to match the best gaussian to avoid duplicates, this makes 
    %the assumption that if two proteins have a very small Eucldian distance the curve fitting give 
    %nearly identical guassian output
for Eucldian_counter2=1: Dimensions_Gaus_import
    if (Interactions_Based_only_Eucldian_distance(Eucldian_counter1,Eucldian_counter2) == 1)
    %Note we have limited the interactions to those with an Eucledian
    %distance which gives a precision over 0.8 and with a center within 2
    %as the guassian measurement defines a single point while eucledian
    %distance is measurement over all points.
    Eu_distance_matched(internal_counter,1)=Protein.Isoform(Eucldian_counter1); %Write out name of protein 1
    Eu_distance_matched{internal_counter,2}=Protein.Center(Eucldian_counter1);
    Eu_distance_matched(internal_counter,3)=Protein.Isoform(Eucldian_counter2); %Write out name of protien 2
    Eu_distance_matched{internal_counter,4}=Protein.Center(Eucldian_counter2);
    %protein 1, determine how gaussians were detected for protein 1
    [L1, L1] = ismember(Protein.Isoform,Protein.Isoform(Eucldian_counter1)); 
    Eu_distance_matched{internal_counter,5}= sum(L1);
    %protein 2, determine how gaussians were detected for protein 2
    [L2, L2] = ismember(Protein.Isoform,Protein.Isoform(Eucldian_counter2)); 
    Eu_distance_matched{internal_counter,6}= sum(L1);
    %Is the interaction in the Corum database
    if  TP_Matrix(Eucldian_counter1,Eucldian_counter2)== 1 %Determine if interaction is known in Corum
        Eu_distance_matched{internal_counter,7}= 1; %if true place 1 in column seven
    else
        Eu_distance_matched{internal_counter,7}= 0; %if false place 0 in column seven
    end
    %Are both proteins in the Corum database
    if  Int_matrix(Eucldian_counter1,Eucldian_counter2)== 1 %Determine both proteins in Corum
        Eu_distance_matched{internal_counter,8}= 1; %if true place 1 in column eight
    else
        Eu_distance_matched{internal_counter,8}= 0; %if false place 0 in column eight
    end
    %Add true (1) to matrix and counter
    Final_Interactions_Based_only_Eucldian_distance(Eucldian_counter1,Eucldian_counter2)= 1;
    internal_counter=internal_counter+1;
    end
end
end

Unique_Interactions_Based_only_Eucldian_distance=triu(Final_Interactions_Based_only_Eucldian_distance);
Total_Interactions_Based_only_Eucldian_distance(replicate_counter,1)= sum(Final_Interactions_Based_only_Eucldian_distance(:));
Count_Interactions_Based_only_Eucldian_distance(replicate_counter,1)= sum(Unique_Interactions_Based_only_Eucldian_distance(:));

%define Eu distance file name
Interactions_Eu_name=strcat('Interactions_determined_based_on_Euclidean_distance_replicate',mat2str(replicate_counter),'.csv');

fid_Eu= fopen(Interactions_Eu_name,'wt'); % Eu dis interaction output 
fprintf (fid_Eu,'%s,%s,%s,%s,%s,%s,%s,%s,\n',... %header for output
                'Protein 1', 'Centre Protein 1', 'Protein 2', 'Centre Protein 2', 'Number of Gaussians Protein 1',...
                'Number of Gaussians Protein 2', 'Interaction Known in Corum', 'Both Proteins in Corum'); %Write Header
for Xiii=1:Total_Interactions_Based_only_Eucldian_distance(replicate_counter,1)   
fprintf(fid_Eu,'%s,%6.3f,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',Eu_distance_matched{Xiii,1}, Eu_distance_matched{Xiii,2},...
    Eu_distance_matched{Xiii,3}, Eu_distance_matched{Xiii,4}, Eu_distance_matched{Xiii,5}, Eu_distance_matched{Xiii,6},...
    Eu_distance_matched{Xiii,7}, Eu_distance_matched{Xiii,8});
end            
fclose(fid_Eu);
cd('../');

%ideally want the widthest Euc, Width, Hight and smallest center that gives
%the most interactions within a defined precision thus we need to determine
%the precision for the range of centers and Euc.

%% Optimise parameters, Determine the number of interaction as center and Euc change and the precision of these interactions
Height_Distance=20.0; % As height has very little effect on interaction determination, this has been set very large
Width_Distangce =20.0; % As width has very little effect on interaction determination, this has been set very large
MatrixHeight = Dist.Height < Height_Distance;
MatrixWidth = Dist.Width < Width_Distangce;

%Determine dimensions of optimisation matrix
%determine the Euc with a FPR of 0.5%
import_name1=strcat('Euc_replicate',mat2str(replicate_counter),'.csv');

%Import results from initial analysis
cd('ROC analysis related Figures'); % move into the ROC figure folder to import Euc file
Eucldian_distance_import = importdata(import_name1);%import Eucldian distance values
cd('../');

%Determine interactions for a FPR of 0.5%
[~, idx1] = min(abs(Eucldian_distance_import.data(:,10) - 0.05)); %determines the postion which gives a FPR of 0.5%
Eucldian_distance_with_0_05FPR= Eucldian_distance_import.data(idx1,1); %determines the Eucldian distance which gives a FPR of 0.5%

%Round Euc distance for a FPR of 0.5%
Eucldian_distance_with_0_05FPR=ceil(Eucldian_distance_with_0_05FPR);

%Define number of data points
Optimisation_Dimension=(Eucldian_distance_with_0_05FPR/0.025);

optimisation_matrix_number_interaction= zeros(Optimisation_Dimension,Optimisation_Dimension);
optimisation_matrix_precision= zeros(Optimisation_Dimension,Optimisation_Dimension);
Euclidian_Distance_optimisation = 0.025:0.025:Eucldian_distance_with_0_05FPR;
Gaussian_fit_Distance_optimisation = 0.025:0.025:Eucldian_distance_with_0_05FPR;

parfor iii=1:Optimisation_Dimension
    for iv=1:Optimisation_Dimension
        NS1=Euclidian_Distance_optimisation(iii);
        NS2=Gaussian_fit_Distance_optimisation(iv);
        MatrixGaussian = Dist.Gaussian_fits < NS1;
        Matrix_Distance = Dist.Euc < NS2;
        Interactions_optimisation_counter = (MatrixGaussian & MatrixHeight & MatrixWidth & Matrix_Distance & inverse_self);
        total_interaction_number= length(find(triu(Interactions_optimisation_counter==1))); % determines the number of interactions in the Upper triangular part of matrix
        optimisation_matrix_number_interaction(iii,iv)=total_interaction_number;
        %Calculate true positives
        TP_Matrix_Int_optimisation = (Interactions_optimisation_counter & TP_Matrix);
        TP_optimisation= length(find(TP_Matrix_Int_optimisation==1));
        %Calculate false positives
        Int2_matrix_optimisation= (Interactions_optimisation_counter & Int_matrix);
        Int3_matrix_optimisation= (Int2_matrix_optimisation & inverse_self);
        FP_optimisation= length(find(Int3_matrix_optimisation==1))-TP_optimisation;
        %Precision
        Precision = TP_optimisation/(TP_optimisation+FP_optimisation);
        optimisation_matrix_precision(iii,iv)=Precision;
    end
end

%write out output of optimisation check
cd('ROC analysis data');
optimisation_matrix_interaction_name=strcat('optimisation_matrix_interaction_numbers_replicate',mat2str(replicate_counter),'.csv');
optimisation_matrix_precision_name=strcat('optimisation_matrix_precision_numbers_replicate',mat2str(replicate_counter),'.csv');
csvwrite(optimisation_matrix_interaction_name,optimisation_matrix_number_interaction);
csvwrite(optimisation_matrix_precision_name,optimisation_matrix_precision);

%% Determine precision final precision across a range of point precision
minimum_precision=(ceil(min(min(optimisation_matrix_precision))*100))/100;
test_precision = minimum_precision:0.01:1;
test_precision =fliplr(test_precision);
number_test_precisions=length(test_precision);
Recall_replicate_level = zeros(1,number_test_precisions);
Precision_replicate_level = zeros(1,number_test_precisions);
TPR_replicate_level = zeros(1,number_test_precisions);
FPR_replicate_level = zeros(1,number_test_precisions);
No_Int_replicate_level= zeros(1,number_test_precisions);

for test_precision_counter=1:number_test_precisions
    
Centre_Euc_require_precision = optimisation_matrix_precision > test_precision(test_precision_counter);   
precision_matrix =optimisation_matrix_number_interaction;
precision_counter= (Optimisation_Dimension*Optimisation_Dimension); %determined by the number of parameters checked during the optimisation 
parfor v=1:precision_counter   % merge the results of checking the precision and number of interactions
        if Centre_Euc_require_precision(v)==1
                precision_matrix(v) = optimisation_matrix_number_interaction(v);
                    else
                       precision_matrix(v)=NaN;
        end             
end

Average_interactions_columns=nanmean(precision_matrix); %Use to measure what the PPI landscape as Center varies with desired precision 
Average_interactions=nanmean(Average_interactions_columns); %Use to measure what the global PPI landscape over Euc and Center

maximum_counter=Optimisation_Dimension; %define based on the optimisation_matrix number 
Interactions_based_one_summed_precisions = zeros(Dimensions_Gaus_import,Dimensions_Gaus_import);
parfor vi=1:maximum_counter   % scan down the Euc row and find which values results the required precision and maximum number of interactions
    Maximum_interactions_at_given_Center = max(precision_matrix(vi,:));
    if ~isnan(Maximum_interactions_at_given_Center)
    inds =find(precision_matrix(vi,:)==Maximum_interactions_at_given_Center);
    [col] = inds(1); %this takes the first number of inds thus ensure the most conservative Center will be used
    optimised_center = Gaussian_fit_Distance_optimisation(vi);
    optimised_Euc = Euclidian_Distance_optimisation(col);
    
    %use optimisation code to maximum out
    MatrixGaussian = Dist.Gaussian_fits < optimised_center;
    Matrix_Distance = Dist.Euc < optimised_Euc;

    % determine final interaction matrices
    Interactions__precisions = (MatrixGaussian & MatrixHeight & MatrixWidth & Matrix_Distance & inverse_self);
    Interactions_based_one_summed_precisions = (Interactions_based_one_summed_precisions | Interactions__precisions);
    end
    
end

% determine final interaction matrices
Final_interactions = (Final_Interactions_Based_only_Eucldian_distance | Interactions_based_one_summed_precisions); %test just interactions on Euc dist

        % Calculate True possitives
        TP_Matrix_Int = (Final_interactions & TP_Matrix);
        TP= length(find(TP_Matrix_Int==1));

        % Calculate false negatives
        MaxTP= length(find(TP_Matrix==1));
        FN= MaxTP-TP;

        % Calculate false positives
        Int2_matrix= (Final_interactions & Int_matrix);
        Int3_matrix= (Int2_matrix & inverse_self);
        FP= length(find(Int3_matrix==1))-TP;
        
        % Find True negatives
        TN= (length(find(Int_matrix==1)))-FP;


%% results of precision analysis within replicate
Recall_replicate_level(test_precision_counter) = TP/(TP+FN);
Precision_replicate_level(test_precision_counter) = TP/(TP+FP);
TPR_replicate_level(test_precision_counter) = TP/(TP+FN);
FPR_replicate_level(test_precision_counter) =FP/(FP+TN);
No_Int_replicate_level(test_precision_counter) = length(find(triu(Final_interactions==1)));

end

Precision_withing_replicate_summary=strcat('Summed_Precision_Summary_replicate',mat2str(replicate_counter),'.csv');
Summed_Precision= fopen(Precision_withing_replicate_summary,'wt'); % create the summary file of the interaction output 
fprintf (Summed_Precision,'%s,%s,%s,%s,%s,%s,\n',... %header for OutputGaus output
               'Precision', 'Summed Recall', 'Summed Precision', 'Summed TPR', 'Summed FPR', 'Summed No_Int'); %Write Header
for Precision_counter=1:number_test_precisions
fprintf (Summed_Precision,'%6.4f, %6.4f,%6.4f,%6.4f,%6.4f,%6.4f,\n',...
                  test_precision(Precision_counter),  Recall_replicate_level(Precision_counter),...
                  Precision_replicate_level(Precision_counter), TPR_replicate_level(Precision_counter),...
                  FPR_replicate_level(Precision_counter), No_Int_replicate_level(Precision_counter));
end
fclose(Summed_Precision);

%Determine interactions within precision of >0.70 based on Eu distance
[~, idx] = min(abs(Precision_replicate_level - 0.7)); %determines the postion which gives a precision closest to 0.7
distance_with_precision_abover_0_7= test_precision(1,idx); %determines the point precision which gives a precision closest to 0.7

%% Define precision required
%note This module has been to designed to use the settings which give the highest number of high precision interactions with the widest Euc
Required_precision = distance_with_precision_abover_0_7; %set precision
Centre_Euc_with_require_precision = optimisation_matrix_precision > Required_precision;
Result_matrix =optimisation_matrix_number_interaction;
precision_counter= Optimisation_Dimension*Optimisation_Dimension; %determined by the number of parameters checked during the optimisation, thus if 200*200 matrix this number is 40000
parfor v=1:precision_counter   % merge the results of checking the precision and number of interactions
        if Centre_Euc_with_require_precision(v)==1
                Result_matrix(v) = optimisation_matrix_number_interaction(v);
                    else
                       Result_matrix(v)=NaN;
        end             
end
Average_interactions_columns=nanmean(Result_matrix); %Use to measure what the PPI landscape as Center varies with desired precision 
Average_interactions=nanmean(Average_interactions_columns); %Use to measure what the global PPI landscape over Euc and Center

optimisation_matrix_required_precision_name=strcat('optimisation_matrix_interaction_numbers_with_required_precision_70pc_replicate',mat2str(replicate_counter),'.csv');
csvwrite(optimisation_matrix_required_precision_name,Result_matrix); %write out the result matrix for use later

%determine the Euc and Center which give the user defined precision

maximum_counter=Optimisation_Dimension; %define based on the optimisation_matrix number 
Interactions_based_one_summed_precisions = zeros(Dimensions_Gaus_import,Dimensions_Gaus_import);
parfor vi=1:maximum_counter   % scan down the Euc row and find which values results the required precision and maximum number of interactions
    Maximum_interactions_at_given_Center = max(Result_matrix(vi,:));
    if ~isnan(Maximum_interactions_at_given_Center)
    inds =find(Result_matrix(vi,:)==Maximum_interactions_at_given_Center);
    [col] = inds(1); %this takes the first number of inds thus ensure the most conservative Center will be used
    optimised_center = Gaussian_fit_Distance_optimisation(vi);
    optimised_Euc = Euclidian_Distance_optimisation(col);
    
    %use optimisation code to maximum out
    MatrixGaussian = Dist.Gaussian_fits < optimised_center;
    Matrix_Distance = Dist.Euc < optimised_Euc;

    % determine final interaction matrices
    Interactions__precisions = (MatrixGaussian & MatrixHeight & MatrixWidth & Matrix_Distance & inverse_self);
    Interactions_based_one_summed_precisions = (Interactions_based_one_summed_precisions | Interactions__precisions);
    end
end

% determine final interaction matrices
Final_interactions = (Final_Interactions_Based_only_Eucldian_distance | Interactions_based_one_summed_precisions); %test just interactions on Euc dist

        % Calculate True possitives
        TP_Matrix_Int = (Final_interactions & TP_Matrix);
        TP= length(find(TP_Matrix_Int==1));

        % Calculate false negatives
        MaxTP= length(find(TP_Matrix==1));
        FN= MaxTP-TP;

        % Calculate false positives
        Int2_matrix= (Final_interactions & Int_matrix);
        Int3_matrix= (Int2_matrix & inverse_self);
        FP= length(find(Int3_matrix==1))-TP;
        
        % Find True negatives
        TN= (length(find(Int_matrix==1)))-FP;
        
        %Generate negative matrix for final interaction analysis
        Negative1_matrix=(Int_matrix  & ~Final_interactions);
        Negative2_matrix=(Negative1_matrix & inverse_self);  

%% final results
Recall = TP/(TP+FN)
Precision = TP/(TP+FP)
TPR = TP/(TP+FN)
FPR =FP/(FP+TN)
No_Int= length(find(triu(Final_interactions==1)))


Summary_Results_70pc_name=strcat('Summary_Results_70pc_replicate',mat2str(replicate_counter),'.csv');
fid_70_2= fopen(Summary_Results_70pc_name,'wt'); % create the summary file of the interaction output 
fprintf (fid_70_2,'%s,%s,%s,%s,%s,\n',... %header for OutputGaus output
                'Recall', 'Precision', 'TPR', 'FPR', 'No_Int'); %Write Header
fprintf (fid_70_2,'%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,\n',...
                    Recall, Precision, TPR, FPR, No_Int);
fprintf (fid_70_2,',\n');             
fclose(fid_70_2);

%% write out final interaction matrix for use to determine the global FDR, TPR, recall and precision across the experiments
Determined_interactions_70pc_precision_name=strcat('Determined_interactions_70pc_precision_replicate',mat2str(replicate_counter),'.csv');
csvwrite(Determined_interactions_70pc_precision_name,Final_interactions); %write out the result matrix for use later

%% Find "Binary" interactions

Protein_elution= strcat(num2str(C),'*',Protein.Isoform);
[n,m]=size(Dist.Euc);
Interactions =cell(1,1);

Protein_interactions=cell(1,2);
Protein_interactions_center=zeros(1,2);
Delta_EucDist=zeros(1,1);
Delta_Center=zeros(1,1);
Delta_Height=zeros(1,1);
Delta_Width=zeros(1,1);
Both_proteins_Corum=zeros(1,1);
Interaction_in_Corum=zeros(1,1);
final_count=0;
U = triu(Final_interactions);    %Take half of matrix therefor unique interactions
for vii=1:n
    for viii= 1:m
        if U(vii,viii)==1
            Int1=Protein_elution(vii);
            Int2=Protein_elution(viii);
            if strcmp(Int1, Int2)
            else
                final_count=1+final_count;
                Interactions{final_count,1}= strcat(Int1,'-',Int2);
                Delta_EucDist(final_count,1)= Dist.Euc(vii,viii);
                Delta_Center(final_count,1)= abs(C(vii)-C(viii));
                Delta_Height(final_count,1)= abs(H(vii)-H(viii));
                Delta_Width(final_count,1)= abs(W(vii)-W(viii));
                Both_proteins_Corum(final_count,1)= Int_matrix(vii,viii);
                Interaction_in_Corum(final_count,1)= TP_Matrix(vii,viii);
                Protein_interactions(final_count,1)=Protein.Isoform(vii);
                Protein_interactions(final_count,2)=Protein.Isoform(viii);
                Protein_interactions_center(final_count,1)=C(vii);
                Protein_interactions_center(final_count,2)=C(viii);
            end
        else
        end   
        
    end   
end

Num_Inter=length(Interactions);
QQ=([Interactions{:}])';

%Print interactions
Output_Interactions_70pc_precision_name=strcat('Output_Interactions_70pc_precision_replicate',mat2str(replicate_counter),'.csv');
fid11 = fopen(Output_Interactions_70pc_precision_name,'w');
fprintf (fid11,'%s,%s,%s,%s,%s,%s,%s,\n',...
        'Interaction', 'Delta Height', 'Delta Center', 'Delta Width', 'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for ix=1:Num_Inter    
fprintf(fid11,'%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',QQ{ix}, Delta_Height(ix), Delta_Center(ix), Delta_Width(ix), Delta_EucDist(ix), Both_proteins_Corum(ix), Interaction_in_Corum(ix));
end
fclose(fid11);

List_Interactions_70pc_precision_name=strcat('Interactions_list_70pc_precision_replicate',mat2str(replicate_counter),'.csv');
fid11B = fopen(List_Interactions_70pc_precision_name,'w');
fprintf (fid11B,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
        'Replicate','Protein A','Protein B', 'Center A', 'Center B',...
        'Delta Height', 'Delta Center', 'Delta Width',...
        'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for ix=1:Num_Inter    
fprintf(fid11B,'%6f,%s,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',replicate_counter,...
    Protein_interactions{ix,1},Protein_interactions{ix,2},Protein_interactions_center(ix,1),Protein_interactions_center(ix,2),...
    Delta_Height(ix),Delta_Center(ix), Delta_Width(ix), Delta_EucDist(ix),...
    Both_proteins_Corum(ix), Interaction_in_Corum(ix));
end
fclose(fid11B);

%% Create list of predicted negative interactions (FN +TN)
%Determine the false "binary" interaction list
U1=triu(Negative2_matrix); %Take half of matrix therefor unique interactions

%change directory
cd('../');
name_negative_file=strcat('Negative interactions_70pc_replicate_',mat2str(replicate_counter));
mkdir(name_negative_file);
parfor vii=1:n
    for viii= 1:m
        if U1(vii,viii)==1
            Int1=Protein_elution(vii);
            Int2=Protein_elution(viii);
            if strcmp(Int1, Int2)
            else
                FN_Interactions= {strcat(Int1,'-',Int2)};
                FN_Delta_EucDist= Dist.Euc(vii,viii);
                FN_Delta_Center= abs(C(vii)-C(viii));
                FN_Delta_Height= abs(H(vii)-H(viii));
                FN_Delta_Width= abs(W(vii)-W(viii));
                FN_Both_proteins_Corum= Int_matrix(vii,viii);
                FN_Interaction_in_Corum= TP_Matrix(vii,viii);
                FN_Protein_interactions1=Protein.Isoform(vii);
                FN_Protein_interactions2=Protein.Isoform(viii);
                FN_Protein_interactions_center1=C(vii);
                FN_Protein_interactions_center2=C(viii);
                
                %Write out data file
                fileName = strcat(mat2str(vii),'_',mat2str(viii),'_Interaction','.csv');
                fid1A = fopen(fileName,'at');
                fprintf(fid1A,'%6f,%s,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',...
                replicate_counter,FN_Protein_interactions1{1},...
                FN_Protein_interactions2{1},FN_Protein_interactions_center1(1),...
                FN_Protein_interactions_center2(1), FN_Delta_Height(1), FN_Delta_Center(1),...
                FN_Delta_Width(1),FN_Delta_EucDist(1), FN_Both_proteins_Corum(1), FN_Interaction_in_Corum(1));
                fclose(fid1A);
                movefile(fileName, name_negative_file);
            end
        else
        end   
        
    end   
end


%read dictory and concatinate negative interactions into a single file
dirData2 = dir(name_negative_file);      %# Get directory contents
dirData2 = dirData2(~[dirData2.isdir]);  %# Use only the file data
negative_fileNames = {dirData2.name};           %# Get file names

ExpName1 = strcat('^(\d+)\_(\d+)\','_Interaction','.csv');
Incorum_Neg_fileData = regexp(negative_fileNames,ExpName1,'tokens');  %# Find tokens
index1 = ~cellfun('isempty',Incorum_Neg_fileData);        %# Find index of matches
negative_fileNames = negative_fileNames(index1);        %# remove empty cells
Incorum_Neg_fileData = [Incorum_Neg_fileData{index1}];                %# Remove non-matches
Incorum_Neg_fileData = vertcat(Incorum_Neg_fileData{:});             %# Format token data
              
%sorting protein interactions in Corum in correct numerical order
nFiles1 = size(Incorum_Neg_fileData,1);              %# Number of files matching format
Incorum_Neg_fileData = str2double(Incorum_Neg_fileData);        %# Convert from strings to numbers
[Incorum_Neg_fileData,index1] = sortrows(Incorum_Neg_fileData,1:2);  %# Sort numeric values
negative_fileNames =  negative_fileNames(index1);               %# Apply sort to file names

Number_of_Neg_interactions = length(negative_fileNames);
Combined_Neg_interactions_data=cell(Number_of_Neg_interactions,11);

for xiv = 1:Number_of_Neg_interactions 
    %load file of interest
    fileName7 = strcat(name_negative_file,'\',negative_fileNames(xiv));
    File_to_be_read = fileName7{1};
    
    %import data files
    f = fopen (File_to_be_read);
    Processed_data = textscan(f, '%s %s %s %s %s %s %s %s %s %s %s %s %s', 'Delimiter', ',');
    fclose(f);  

    %Copy values to arrays to use in downstream analysis
    %Value 1, As the first value is the replicate there is not need to save
    %value
    
    %Value 2
    prot1_protein_names_entry=Processed_data{2};
    Neg_Corum_prot1_protein_names{xiv}=prot1_protein_names_entry{1};
    %Value 3
    prot2_protein_names_entry=Processed_data{3};
    Neg_Corum_prot2_protein_names{xiv}=prot2_protein_names_entry{1};
    %Value 4
    prot1_center=Processed_data{4};
    Neg_Corum_prot1_center(xiv)=str2double(prot1_center(1));
    %Value 5
    prot2_center=Processed_data{5};
    Neg_Corum_prot2_center(xiv)=str2double(prot2_center{1});
    %Value 6
    Neg_Delta_height_value=Processed_data{6};
    Neg_Delta_height(xiv)=str2double(Neg_Delta_height_value{1});
    %Value 7
    Neg_Center_value=Processed_data{7};
    Neg_Center(xiv)=str2double(Neg_Center_value(1));
    %Value 8
    Neg_Width_value=Processed_data{8};
    Neg_Width(xiv)=str2double(Neg_Width_value(1));
    %Value 9
    Neg_Euc_value=Processed_data{9};
    Neg_Euc(xiv)=str2double(Neg_Euc_value{1});
    %Value 10
    Neg_both_in_corum_value=Processed_data{10};
    Neg_both_in_corum(xiv)=str2double(Neg_both_in_corum_value{1});
    %Value 11
    Neg_interaction_in_corum_value=Processed_data{11};
    Neg_interaction_in_corum(xiv)=str2double(Neg_interaction_in_corum_value(1));
end


%Write out negative interaction list, used to create the final TPR, FPR
%precision and recall values
cd('ROC analysis data');
Neg_proteins_interactions_name=strcat('Neg_proteins_interactions_70pc_replicate_',mat2str(replicate_counter),'.csv');
fid_Neg= fopen(Neg_proteins_interactions_name,'w'); 
fprintf (fid_Neg,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
        'Replicate','Protein A','Protein B', 'Center A', 'Center B',...
        'Delta Height', 'Delta Center', 'Delta Width',...
        'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for xv=1:Number_of_Neg_interactions   
fprintf(fid_Neg,'%6f,%s,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',...
replicate_counter(1) ,Neg_Corum_prot1_protein_names{xv}, Neg_Corum_prot2_protein_names{xv},...
Neg_Corum_prot1_center(xv),Neg_Corum_prot2_center(xv), Neg_Delta_height(xv),...
Neg_Center(xv),Neg_Width(xv),Neg_Euc(xv),Neg_both_in_corum(xv),Neg_interaction_in_corum(xv));
end     
fclose(fid_Neg);

%% List of interactions
List_of_final_interactions = Final_interactions;    %As we want the complete list do not use triu
Protein_elution_1= strcat(num2str(C),'*',Protein.Isoform);
[n,m]=size(Dist.Euc);
row_counter=0;
empty_cell=cell(1,1);
for vii=1:n
    List_of_final_interactions_Protein.Guassian(vii,1)= Protein_elution_1(vii);
    List_of_final_interactions_Protein.Protein(vii,1)= Protein.Isoform(vii);
    List_of_final_interactions_Protein.Protein_NoIsoform(vii,1)= Protein.NoIsoform(vii);
        for MajorID_protein_counter = 1:dimension_Protein_MajorID_NoIsoforms(2) 
        List_of_final_interactions_Protein.Protein_MajorID_NoIsoforms(vii,MajorID_protein_counter)= Protein.MajorID_NoIsoforms(vii,MajorID_protein_counter);
        end
    List_of_final_interactions_Protein.Center(vii,1)= C(vii);
    final_count=0;
    for viii= 1:m
        if List_of_final_interactions(vii,viii)== 1 %Determine if proteins interact
            final_count=1+final_count;
            List_of_final_interactions_Protein.Interaction_Gaussian(vii,final_count)=Protein_elution_1(viii);
            List_of_final_interactions_Protein.Interaction_Protein(vii,final_count)=Protein.Isoform(viii);
            List_of_final_interactions_Protein.Interaction_Protein_NoIsoform(vii,final_count)=Protein.NoIsoform(viii);
            List_of_final_interactions_Protein.Interaction_Center(vii,final_count)=C(viii);
            if  TP_Matrix(vii,viii)== 1 %Determine if interaction is known in Corum
                List_of_final_interactions_Protein.Interaction_known_in_corum(vii,final_count)= 1;
            else
                List_of_final_interactions_Protein.Interaction_known_in_corum(vii,final_count)= 0;
            end
        end
    end
    if final_count == 0
    List_of_final_interactions_Protein.Interaction_Gaussian(vii,1)= {NaN};
    List_of_final_interactions_Protein.Interaction_Protein(vii,1)={NaN};
    List_of_final_interactions_Protein.Interaction_Center(vii,1)= 0;
    List_of_final_interactions_Protein.Interaction_known_in_corum(vii,1)=0;
    end
end

%Remove values in List of final with NaN
List_of_final_interactions_Protein.Interaction_Gaussian(cellfun(@(x) any(isnan(x)),List_of_final_interactions_Protein.Interaction_Gaussian)) = {''};
List_of_final_interactions_Protein.Interaction_Protein(cellfun(@(x) any(isnan(x)),List_of_final_interactions_Protein.Interaction_Protein)) = {''};

%Determine the dimension of List_of_final_interactions_Protein.Interaction_Protein array
Maximum_observed_interactions=size(List_of_final_interactions_Protein.Interaction_Protein); % determine the dimension of the interaction array in structure

%Count number of corum interactions and total interactions
List_of_final_interactions_Protein.List_Interaction_Gaussians=cell(Maximum_observed_interactions(1),1);
List_of_final_interactions_Protein.Total_Interaction=cell(Maximum_observed_interactions(1),1);
for Interaction_number_counter1=1:Maximum_observed_interactions(1)
    Interaction_counter_for_protein=0;
    for Interaction_number_counter2=1:Maximum_observed_interactions(2)
    Guassian_to_combine(Interaction_number_counter2)=List_of_final_interactions_Protein.Interaction_Gaussian(Interaction_number_counter1,Interaction_number_counter2);
    if any(Guassian_to_combine{Interaction_number_counter2})== 1 % test if value is not zero
        String_name= strcat(Guassian_to_combine{Interaction_number_counter2},' ; ');
        List_of_final_interactions_Protein.List_Interaction_Gaussians{Interaction_number_counter1}=strcat(String_name,List_of_final_interactions_Protein.List_Interaction_Gaussians{Interaction_number_counter1});
        Interaction_counter_for_protein=Interaction_counter_for_protein+1;
    end
    end
    List_of_final_interactions_Protein.Total_Interaction{Interaction_number_counter1,1}=Interaction_counter_for_protein;
end

%Sum number of known interactions in corum for protein
for Interaction_number_counter1=1:Maximum_observed_interactions(1)
for Interaction_number_counter2=1:Maximum_observed_interactions(2)
Detected_Center_Interaction_known_in_corum(Interaction_number_counter2)= List_of_final_interactions_Protein.Interaction_known_in_corum(Interaction_number_counter1,Interaction_number_counter2);
end
List_of_final_interactions_Protein.Number_Interaction_known_in_corum(Interaction_number_counter1,1)= sum(Detected_Center_Interaction_known_in_corum);
end

%Protein list of interactions
List_of_final_interactions_Protein.List_Interaction_Protein=cell(Maximum_observed_interactions(1),1);
for protein_table1=1:Maximum_observed_interactions(1)
    for protein_table2=1:Maximum_observed_interactions(2)
    Proteins_to_combine(protein_table2)=List_of_final_interactions_Protein.Interaction_Protein(protein_table1,protein_table2);
    if any(Proteins_to_combine{protein_table2})== 1 % test if value is not zero
        String_name= strcat(Proteins_to_combine{protein_table2},' ; ');
        List_of_final_interactions_Protein.List_Interaction_Protein{protein_table1,1}=strcat(String_name,List_of_final_interactions_Protein.List_Interaction_Protein{protein_table1});
    end
    end
end

%Copy List_of_final_interactions_Protein_70pc to use for other experiuments
List_of_final_interactions_Protein_70pc=List_of_final_interactions_Protein;

%Sort structure by Center/Gaussian write out table
[location2,location2]=sort([List_of_final_interactions_Protein.Center]); %Sort strcuture based on protein name
for sort_counter1=1:Maximum_observed_interactions(1)
Sorted_of_final_interactions_Center_70pc.Guassian{sort_counter1,1} = List_of_final_interactions_Protein.Guassian{location2(sort_counter1)};
Sorted_of_final_interactions_Center_70pc.Protein{sort_counter1,1}= List_of_final_interactions_Protein.Protein{location2(sort_counter1)};
Sorted_of_final_interactions_Center_70pc.Protein_NoIsoform{sort_counter1,1}=List_of_final_interactions_Protein.Protein_NoIsoform{location2(sort_counter1)};
    for MajorID_protein_counter = 1:dimension_Protein_MajorID_NoIsoforms(2)
    Sorted_of_final_interactions_Center_70pc.Protein_MajorID_NoIsoforms{sort_counter1,MajorID_protein_counter}=List_of_final_interactions_Protein.Protein_MajorID_NoIsoforms{location2(sort_counter1),MajorID_protein_counter};    
    end
Sorted_of_final_interactions_Center_70pc.Center{sort_counter1,1}= List_of_final_interactions_Protein.Center(location2(sort_counter1));
    for sort_counter2=1:Maximum_observed_interactions(2)
    Sorted_of_final_interactions_Center_70pc.Interaction_Gaussian{sort_counter1,sort_counter2}= List_of_final_interactions_Protein.Interaction_Gaussian{location2(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Center_70pc.Interaction_Protein{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Protein{location2(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Center_70pc.Interaction_Center{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Center(location2(sort_counter1),sort_counter2);
    Sorted_of_final_interactions_Center_70pc.Interaction_known_in_corum{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_known_in_corum(location2(sort_counter1),sort_counter2);
    end
Sorted_of_final_interactions_Center_70pc.Number_Interaction_known_in_corum(sort_counter1,1)= List_of_final_interactions_Protein.Number_Interaction_known_in_corum(location2(sort_counter1));
Sorted_of_final_interactions_Center_70pc.Total_Interaction(sort_counter1,1)=List_of_final_interactions_Protein.Total_Interaction(location2(sort_counter1));
Sorted_of_final_interactions_Center_70pc.List_Interaction_Gaussians(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Gaussians(location2(sort_counter1));
Sorted_of_final_interactions_Center_70pc.List_Interaction_Protein(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Protein(location2(sort_counter1));
end

%write out sorted structure array, Gaussian table
Sort_Interactions_70pc_precision_name=strcat('Sorted_Output_Interactions_Gaussians_70pc_replicate',mat2str(replicate_counter),'.csv');
fid_sorted_interactions_70_1 = fopen(Sort_Interactions_70pc_precision_name,'w');
fprintf (fid_sorted_interactions_70_1,'%s,%s,%s,%s,%s,%s,%s,\n',...
'Guassian', 'Protein_name(Isoform)','Protein_name(NoIsoform)','Center','Number of interactions','Number of interactions in Corum','Protein Interactions determined');               %Write Header
for write_out_sorted=1:Maximum_observed_interactions(1)
       fprintf(fid_sorted_interactions_70_1, '%s,%s,%s,%6.4f,%6.4f,%6.4f,%s,\n',...
       Sorted_of_final_interactions_Center_70pc.Guassian{write_out_sorted},...
       Sorted_of_final_interactions_Center_70pc.Protein{write_out_sorted},...
       Sorted_of_final_interactions_Center_70pc.Protein_NoIsoform{write_out_sorted},...
       Sorted_of_final_interactions_Center_70pc.Center{write_out_sorted},...
       Sorted_of_final_interactions_Center_70pc.Total_Interaction{write_out_sorted},...
       Sorted_of_final_interactions_Center_70pc.Number_Interaction_known_in_corum(write_out_sorted),...
       Sorted_of_final_interactions_Center_70pc.List_Interaction_Gaussians{write_out_sorted});
end
fclose(fid_sorted_interactions_70_1);

%Sort strcuture by Protein name
[location,location]=sort([List_of_final_interactions_Protein.Protein]); %Sort strcuture based on protein name
for sort_counter1=1:Maximum_observed_interactions(1)
Sorted_of_final_interactions_Protein_70pc.Guassian{sort_counter1,1} = List_of_final_interactions_Protein.Guassian{location(sort_counter1)};
Sorted_of_final_interactions_Protein_70pc.Protein{sort_counter1,1}= List_of_final_interactions_Protein.Protein{location(sort_counter1)};
Sorted_of_final_interactions_Protein_70pc.Protein_NoIsoform{sort_counter1,1}=List_of_final_interactions_Protein.Protein_NoIsoform{location(sort_counter1)};
    for MajorID_protein_counter = 1:dimension_Protein_MajorID_NoIsoforms(2)
    Sorted_of_final_interactions_Protein_70pc.Protein_MajorID_NoIsoforms{sort_counter1,MajorID_protein_counter}=List_of_final_interactions_Protein.Protein_MajorID_NoIsoforms{location(sort_counter1),MajorID_protein_counter};    
    end
Sorted_of_final_interactions_Protein_70pc.Center{sort_counter1,1}= List_of_final_interactions_Protein.Center(location(sort_counter1));
    for sort_counter2=1:Maximum_observed_interactions(2)
    Sorted_of_final_interactions_Protein_70pc.Interaction_Gaussian{sort_counter1,sort_counter2}= List_of_final_interactions_Protein.Interaction_Gaussian{location(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Protein_70pc.Interaction_Protein{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Protein{location(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Protein_70pc.Interaction_Center{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Center(location(sort_counter1),sort_counter2);
    Sorted_of_final_interactions_Protein_70pc.Interaction_known_in_corum{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_known_in_corum(location(sort_counter1),sort_counter2);
    end
Sorted_of_final_interactions_Protein_70pc.Number_Interaction_known_in_corum(sort_counter1,1)= List_of_final_interactions_Protein.Number_Interaction_known_in_corum(location(sort_counter1));
Sorted_of_final_interactions_Protein_70pc.Total_Interaction(sort_counter1,1)=List_of_final_interactions_Protein.Total_Interaction(location(sort_counter1));
Sorted_of_final_interactions_Protein_70pc.List_Interaction_Gaussians(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Gaussians(location(sort_counter1));
Sorted_of_final_interactions_Protein_70pc.List_Interaction_Protein(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Protein(location(sort_counter1));    
end

%write out sorted structure array, protein table
Sort_Interactions_70pc_Proteins_precision_name=strcat('Sorted_Output_Interactions_Proteins_70pc_replicate',mat2str(replicate_counter),'.csv');
fid_sorted_interactions_70_2 = fopen(Sort_Interactions_70pc_Proteins_precision_name,'w');
fprintf (fid_sorted_interactions_70_2,'%s,%s,%s,%s,%s,%s,%s,\n',...
'Guassian', 'Protein_name (Isoform)','Protein_name (NoIsoform)','Center','Number of interactions','Number of interactions in Corum','Protein Interactions determined');               %Write Header
for write_out_sorted=1:Maximum_observed_interactions(1)
       fprintf(fid_sorted_interactions_70_2, '%s,%s,%s,%6.4f,%6.4f,%6.4f,%s,\n',...
       Sorted_of_final_interactions_Protein_70pc.Guassian{write_out_sorted},...
       Sorted_of_final_interactions_Protein_70pc.Protein{write_out_sorted},...
       Sorted_of_final_interactions_Protein_70pc.Protein_NoIsoform{write_out_sorted},...
       Sorted_of_final_interactions_Protein_70pc.Center{write_out_sorted},...
       Sorted_of_final_interactions_Protein_70pc.Total_Interaction{write_out_sorted},...
       Sorted_of_final_interactions_Protein_70pc.Number_Interaction_known_in_corum(write_out_sorted),...
       Sorted_of_final_interactions_Protein_70pc.List_Interaction_Protein{write_out_sorted});
end
fclose(fid_sorted_interactions_70_2);

%% Define precision required 60pc
[~, idx] = min(abs(Precision_replicate_level - 0.6)); %determines the postion which gives a precision closest to 0.6
distance_with_precision_abover_0_6= test_precision(1,idx); %determines the point precision which gives a precision closest to 0.6

%note This module has been to designed to use the settings which give the highest number of high precision interactions with the widest Euc
Required_precision = distance_with_precision_abover_0_6; %set precision
Centre_Euc_with_require_precision = optimisation_matrix_precision > Required_precision;
Result_matrix =optimisation_matrix_number_interaction;
precision_counter= Optimisation_Dimension*Optimisation_Dimension; %determined by the number of parameters checked during the optimisation, thus if 200*200 matrix this number is 40000
parfor v=1:precision_counter   % merge the results of checking the precision and number of interactions
        if Centre_Euc_with_require_precision(v)==1
                Result_matrix(v) = optimisation_matrix_number_interaction(v);
                    else
                       Result_matrix(v)=NaN;
        end             
end
Average_interactions_columns=nanmean(Result_matrix); %Use to measure what the PPI landscape as Center varies with desired precision 
Average_interactions=nanmean(Average_interactions_columns); %Use to measure what the global PPI landscape over Euc and Center

optimisation_matrix_required_precision_name=strcat('optimisation_matrix_interaction_numbers_with_required_precision_60pc_replicate',mat2str(replicate_counter),'.csv');
csvwrite(optimisation_matrix_required_precision_name,Result_matrix); %write out the result matrix for use later

%determine the Euc and Center which give the user defined precision

maximum_counter=Optimisation_Dimension; %define based on the optimisation_matrix number 
Interactions_based_one_summed_precisions = zeros(Dimensions_Gaus_import,Dimensions_Gaus_import);
parfor vi=1:maximum_counter   % scan down the Euc row and find which values results the required precision and maximum number of interactions
    Maximum_interactions_at_given_Center = max(Result_matrix(vi,:));
    if ~isnan(Maximum_interactions_at_given_Center)
    inds =find(Result_matrix(vi,:)==Maximum_interactions_at_given_Center);
    [col] = inds(1); %this takes the first number of inds thus ensure the most conservative Center will be used
    optimised_center = Gaussian_fit_Distance_optimisation(vi);
    optimised_Euc = Euclidian_Distance_optimisation(col);
    
    %use optimisation code to maximum out
    MatrixGaussian = Dist.Gaussian_fits < optimised_center;
    Matrix_Distance = Dist.Euc < optimised_Euc;

    % determine final interaction matrices
    Interactions__precisions = (MatrixGaussian & MatrixHeight & MatrixWidth & Matrix_Distance & inverse_self);
    Interactions_based_one_summed_precisions = (Interactions_based_one_summed_precisions | Interactions__precisions);
    end
end

% determine final interaction matrices
Final_interactions = (Final_Interactions_Based_only_Eucldian_distance | Interactions_based_one_summed_precisions); %test just interactions on Euc dist

        % Calculate True possitives
        TP_Matrix_Int = (Final_interactions & TP_Matrix);
        TP= length(find(TP_Matrix_Int==1));

        % Calculate false negatives
        MaxTP= length(find(TP_Matrix==1));
        FN= MaxTP-TP;

        % Calculate false positives
        Int2_matrix= (Final_interactions & Int_matrix);
        Int3_matrix= (Int2_matrix & inverse_self);
        FP= length(find(Int3_matrix==1))-TP;
        
        % Find True negatives
        TN= (length(find(Int_matrix==1)))-FP;
        
        %Generate negative matrix for final interaction analysis
        Negative1_matrix=(Int_matrix  & ~Final_interactions);
        Negative2_matrix=(Negative1_matrix & inverse_self); 


%% final results
Recall = TP/(TP+FN)
Precision = TP/(TP+FP)
TPR = TP/(TP+FN)
FPR =FP/(FP+TN)
No_Int= length(find(triu(Final_interactions==1)))

Summary_Results_60pc_name=strcat('Summary_Results_60pc_replicate',mat2str(replicate_counter),'.csv');
fid_60_2= fopen(Summary_Results_60pc_name,'wt'); % create the summary file of the interaction output 
fprintf (fid_60_2,'%s,%s,%s,%s,%s,\n',... %header for OutputGaus output
                'Recall', 'Precision', 'TPR', 'FPR', 'No_Int'); %Write Header
fprintf (fid_60_2,'%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,\n',...
                    Recall, Precision, TPR, FPR, No_Int);
fprintf (fid_60_2,',\n');             
fclose(fid_60_2);

%% write out final interaction matrix for use to determine the global FDR, TPR, recall and precision across the experiments
Determined_interactions_60pc_precision_name=strcat('Determined_interactions_60pc_precision_replicate',mat2str(replicate_counter),'.csv');
csvwrite(Determined_interactions_60pc_precision_name,Final_interactions); %write out the result matrix for use later

%% Find "Binary" interactions

Protein_elution= strcat(num2str(C),'*',Protein.Isoform);
[n,m]=size(Dist.Euc);
Interactions =cell(1,1);

Protein_interactions=cell(1,2);
Protein_interactions_center=zeros(1,2);
Delta_EucDist=zeros(1,1);
Delta_Center=zeros(1,1);
Delta_Height=zeros(1,1);
Delta_Width=zeros(1,1);
Both_proteins_Corum=zeros(1,1);
Interaction_in_Corum=zeros(1,1);
final_count=0;
U = triu(Final_interactions);    %Take half of matrix therefor unique interactions
for vii=1:n
    for viii= 1:m
        if U(vii,viii)==1
            Int1=Protein_elution(vii);
            Int2=Protein_elution(viii);
            if strcmp(Int1, Int2)
            else
                final_count=1+final_count;
                Interactions{final_count,1}= strcat(Int1,'-',Int2);
                Delta_EucDist(final_count,1)= Dist.Euc(vii,viii);
                Delta_Center(final_count,1)= abs(C(vii)-C(viii));
                Delta_Height(final_count,1)= abs(H(vii)-H(viii));
                Delta_Width(final_count,1)= abs(W(vii)-W(viii));
                Both_proteins_Corum(final_count,1)= Int_matrix(vii,viii);
                Interaction_in_Corum(final_count,1)= TP_Matrix(vii,viii);
                Protein_interactions(final_count,1)=Protein.Isoform(vii);
                Protein_interactions(final_count,2)=Protein.Isoform(viii);
                Protein_interactions_center(final_count,1)=C(vii);
                Protein_interactions_center(final_count,2)=C(viii);
            end
        else
        end   
        
    end   
end


Num_Inter=length(Interactions);
QQ=([Interactions{:}])';

%Print interactions
Output_Interactions_60pc_precision_name=strcat('Output_Interactions_60pc_precision_replicate',mat2str(replicate_counter),'.csv');
fid11 = fopen(Output_Interactions_60pc_precision_name,'w');
fprintf (fid11,'%s,%s,%s,%s,%s,%s,%s,\n',...
        'Interaction', 'Delta Height', 'Delta Center', 'Delta Width', 'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for ix=1:Num_Inter    
fprintf(fid11,'%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',QQ{ix}, Delta_Height(ix), Delta_Center(ix), Delta_Width(ix), Delta_EucDist(ix), Both_proteins_Corum(ix), Interaction_in_Corum(ix));
end
fclose(fid11);

List_Interactions_60pc_precision_name=strcat('Interactions_list_60pc_precision_replicate',mat2str(replicate_counter),'.csv');
fid11B = fopen(List_Interactions_60pc_precision_name,'w');
fprintf (fid11B,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
        'Replicate','Protein A','Protein B', 'Center A', 'Center B',...
        'Delta Height', 'Delta Center', 'Delta Width',...
        'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for ix=1:Num_Inter    
fprintf(fid11B,'%6f,%s,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',replicate_counter,... 
    Protein_interactions{ix,1},Protein_interactions{ix,2},Protein_interactions_center(ix,1),Protein_interactions_center(ix,2),...
    Delta_Height(ix),Delta_Center(ix), Delta_Width(ix), Delta_EucDist(ix),...
    Both_proteins_Corum(ix), Interaction_in_Corum(ix));
end
fclose(fid11B);

%% Create list of predicted negative interactions (FN +TN)
%Determine the false "binary" interaction list
U1=triu(Negative2_matrix); %Take half of matrix therefor unique interactions

%change directory
cd('../');
name_negative_file=strcat('Negative interactions_60pc_replicate_',mat2str(replicate_counter));
mkdir(name_negative_file);
parfor vii=1:n
    for viii= 1:m
        if U1(vii,viii)==1
            Int1=Protein_elution(vii);
            Int2=Protein_elution(viii);
            if strcmp(Int1, Int2)
            else
                FN_Interactions= {strcat(Int1,'-',Int2)};
                FN_Delta_EucDist= Dist.Euc(vii,viii);
                FN_Delta_Center= abs(C(vii)-C(viii));
                FN_Delta_Height= abs(H(vii)-H(viii));
                FN_Delta_Width= abs(W(vii)-W(viii));
                FN_Both_proteins_Corum= Int_matrix(vii,viii);
                FN_Interaction_in_Corum= TP_Matrix(vii,viii);
                FN_Protein_interactions1=Protein.Isoform(vii);
                FN_Protein_interactions2=Protein.Isoform(viii);
                FN_Protein_interactions_center1=C(vii);
                FN_Protein_interactions_center2=C(viii);
                
                %Write out data file
                fileName = strcat(mat2str(vii),'_',mat2str(viii),'_Interaction','.csv');
                fid1A = fopen(fileName,'at');
                fprintf(fid1A,'%6f,%s,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',...
                replicate_counter,FN_Protein_interactions1{1},...
                FN_Protein_interactions2{1},FN_Protein_interactions_center1(1),...
                FN_Protein_interactions_center2(1), FN_Delta_Height(1), FN_Delta_Center(1),...
                FN_Delta_Width(1),FN_Delta_EucDist(1), FN_Both_proteins_Corum(1), FN_Interaction_in_Corum(1));
                fclose(fid1A);
                movefile(fileName, name_negative_file);
            end
        else
        end   
        
    end   
end


%read dictory and concatinate negative interactions into a single file
dirData2 = dir(name_negative_file);      %# Get directory contents
dirData2 = dirData2(~[dirData2.isdir]);  %# Use only the file data
negative_fileNames = {dirData2.name};           %# Get file names

ExpName1 = strcat('^(\d+)\_(\d+)\','_Interaction','.csv');
Incorum_Neg_fileData = regexp(negative_fileNames,ExpName1,'tokens');  %# Find tokens
index1 = ~cellfun('isempty',Incorum_Neg_fileData);        %# Find index of matches
negative_fileNames = negative_fileNames(index1);        %# remove empty cells
Incorum_Neg_fileData = [Incorum_Neg_fileData{index1}];                %# Remove non-matches
Incorum_Neg_fileData = vertcat(Incorum_Neg_fileData{:});             %# Format token data
              
%sorting protein interactions in Corum in correct numerical order
nFiles1 = size(Incorum_Neg_fileData,1);              %# Number of files matching format
Incorum_Neg_fileData = str2double(Incorum_Neg_fileData);        %# Convert from strings to numbers
[Incorum_Neg_fileData,index1] = sortrows(Incorum_Neg_fileData,1:2);  %# Sort numeric values
negative_fileNames =  negative_fileNames(index1);               %# Apply sort to file names

Number_of_Neg_interactions = length(negative_fileNames);
Combined_Neg_interactions_data=cell(Number_of_Neg_interactions,11);

for xiv = 1:Number_of_Neg_interactions 
    %load file of interest
    fileName7 = strcat(name_negative_file,'\',negative_fileNames(xiv));
    File_to_be_read = fileName7{1};
    
    %import data files
    f = fopen (File_to_be_read);
    Processed_data = textscan(f, '%s %s %s %s %s %s %s %s %s %s %s %s %s', 'Delimiter', ',');
    fclose(f);  

    %Copy values to arrays to use in downstream analysis
    %Value 1, As the first value is the replicate there is not need to save
    %value
    
    %Value 2
    prot1_protein_names_entry=Processed_data{2};
    Neg_Corum_prot1_protein_names{xiv}=prot1_protein_names_entry{1};
    %Value 3
    prot2_protein_names_entry=Processed_data{3};
    Neg_Corum_prot2_protein_names{xiv}=prot2_protein_names_entry{1};
    %Value 4
    prot1_center=Processed_data{4};
    Neg_Corum_prot1_center(xiv)=str2double(prot1_center(1));
    %Value 5
    prot2_center=Processed_data{5};
    Neg_Corum_prot2_center(xiv)=str2double(prot2_center{1});
    %Value 6
    Neg_Delta_height_value=Processed_data{6};
    Neg_Delta_height(xiv)=str2double(Neg_Delta_height_value{1});
    %Value 7
    Neg_Center_value=Processed_data{7};
    Neg_Center(xiv)=str2double(Neg_Center_value(1));
    %Value 8
    Neg_Width_value=Processed_data{8};
    Neg_Width(xiv)=str2double(Neg_Width_value(1));
    %Value 9
    Neg_Euc_value=Processed_data{9};
    Neg_Euc(xiv)=str2double(Neg_Euc_value{1});
    %Value 10
    Neg_both_in_corum_value=Processed_data{10};
    Neg_both_in_corum(xiv)=str2double(Neg_both_in_corum_value{1});
    %Value 11
    Neg_interaction_in_corum_value=Processed_data{11};
    Neg_interaction_in_corum(xiv)=str2double(Neg_interaction_in_corum_value(1));
end

%Write out negative interaction list, used to create the final TPR, FPR
%precision and recall values
cd('ROC analysis data');
Neg_proteins_interactions_name=strcat('Neg_proteins_interactions_60pc_replicate_',mat2str(replicate_counter),'.csv');
fid_Neg= fopen(Neg_proteins_interactions_name,'w'); 
fprintf (fid_Neg,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
        'Replicate','Protein A','Protein B', 'Center A', 'Center B',...
        'Delta Height', 'Delta Center', 'Delta Width',...
        'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for xv=1:Number_of_Neg_interactions   
fprintf(fid_Neg,'%6f,%s,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',...
replicate_counter(1) ,Neg_Corum_prot1_protein_names{xv}, Neg_Corum_prot2_protein_names{xv},...
Neg_Corum_prot1_center(xv),Neg_Corum_prot2_center(xv), Neg_Delta_height(xv),...
Neg_Center(xv),Neg_Width(xv),Neg_Euc(xv),Neg_both_in_corum(xv),Neg_interaction_in_corum(xv));
end     
fclose(fid_Neg);

%% List of interactions
List_of_final_interactions = Final_interactions;    %As we want the complete list do not use triu
Protein_elution_1= strcat(num2str(C),'*',Protein.Isoform);
[n,m]=size(Dist.Euc);
row_counter=0;
empty_cell=cell(1,1);
for vii=1:n
    List_of_final_interactions_Protein.Guassian(vii,1)= Protein_elution_1(vii);
    List_of_final_interactions_Protein.Protein(vii,1)= Protein.Isoform(vii);
    List_of_final_interactions_Protein.Protein_NoIsoform(vii,1)= Protein.NoIsoform(vii);
        for MajorID_protein_counter = 1:dimension_Protein_MajorID_NoIsoforms(2) 
        List_of_final_interactions_Protein.Protein_MajorID_NoIsoforms(vii,MajorID_protein_counter)= Protein.MajorID_NoIsoforms(vii,MajorID_protein_counter);
        end
    List_of_final_interactions_Protein.Center(vii,1)= C(vii);
    final_count=0;
    for viii= 1:m
        if List_of_final_interactions(vii,viii)== 1 %Determine if proteins interact
            final_count=1+final_count;
            List_of_final_interactions_Protein.Interaction_Gaussian(vii,final_count)=Protein_elution_1(viii);
            List_of_final_interactions_Protein.Interaction_Protein(vii,final_count)=Protein.Isoform(viii);
            List_of_final_interactions_Protein.Interaction_Protein_NoIsoform(vii,final_count)=Protein.NoIsoform(viii);
            List_of_final_interactions_Protein.Interaction_Center(vii,final_count)=C(viii);
            if  TP_Matrix(vii,viii)== 1 %Determine if interaction is known in Corum
                List_of_final_interactions_Protein.Interaction_known_in_corum(vii,final_count)= 1;
            else
                List_of_final_interactions_Protein.Interaction_known_in_corum(vii,final_count)= 0;
            end
        end
    end
    if final_count == 0
    List_of_final_interactions_Protein.Interaction_Gaussian(vii,1)= {NaN};
    List_of_final_interactions_Protein.Interaction_Protein(vii,1)={NaN};
    List_of_final_interactions_Protein.Interaction_Center(vii,1)= 0;
    List_of_final_interactions_Protein.Interaction_known_in_corum(vii,1)=0;
    end
end

%Remove values in List of final with NaN
List_of_final_interactions_Protein.Interaction_Gaussian(cellfun(@(x) any(isnan(x)),List_of_final_interactions_Protein.Interaction_Gaussian)) = {''};
List_of_final_interactions_Protein.Interaction_Protein(cellfun(@(x) any(isnan(x)),List_of_final_interactions_Protein.Interaction_Protein)) = {''};

%Determine the dimension of List_of_final_interactions_Protein.Interaction_Protein array
Maximum_observed_interactions=size(List_of_final_interactions_Protein.Interaction_Protein); % determine the dimension of the interaction array in structure

%Count number of corum interactions and total interactions
List_of_final_interactions_Protein.List_Interaction_Gaussians=cell(Maximum_observed_interactions(1),1);
List_of_final_interactions_Protein.Total_Interaction=cell(Maximum_observed_interactions(1),1);
for Interaction_number_counter1=1:Maximum_observed_interactions(1)
    Interaction_counter_for_protein=0;
    for Interaction_number_counter2=1:Maximum_observed_interactions(2)
    Guassian_to_combine(Interaction_number_counter2)=List_of_final_interactions_Protein.Interaction_Gaussian(Interaction_number_counter1,Interaction_number_counter2);
    if any(Guassian_to_combine{Interaction_number_counter2})== 1 % test if value is not zero
        String_name= strcat(Guassian_to_combine{Interaction_number_counter2},' ; ');
        List_of_final_interactions_Protein.List_Interaction_Gaussians{Interaction_number_counter1}=strcat(String_name,List_of_final_interactions_Protein.List_Interaction_Gaussians{Interaction_number_counter1});
        Interaction_counter_for_protein=Interaction_counter_for_protein+1;
    end
    end
    List_of_final_interactions_Protein.Total_Interaction{Interaction_number_counter1,1}=Interaction_counter_for_protein;
end

%Sum number of known interactions in corum for protein
for Interaction_number_counter1=1:Maximum_observed_interactions(1)
for Interaction_number_counter2=1:Maximum_observed_interactions(2)
Detected_Center_Interaction_known_in_corum(Interaction_number_counter2)= List_of_final_interactions_Protein.Interaction_known_in_corum(Interaction_number_counter1,Interaction_number_counter2);
end
List_of_final_interactions_Protein.Number_Interaction_known_in_corum(Interaction_number_counter1,1)= sum(Detected_Center_Interaction_known_in_corum);
end

%Protein list of interactions
List_of_final_interactions_Protein.List_Interaction_Protein=cell(Maximum_observed_interactions(1),1);
for protein_table1=1:Maximum_observed_interactions(1)
    for protein_table2=1:Maximum_observed_interactions(2)
    Proteins_to_combine(protein_table2)=List_of_final_interactions_Protein.Interaction_Protein(protein_table1,protein_table2);
    if any(Proteins_to_combine{protein_table2})== 1 % test if value is not zero
        String_name= strcat(Proteins_to_combine{protein_table2},' ; ');
        List_of_final_interactions_Protein.List_Interaction_Protein{protein_table1,1}=strcat(String_name,List_of_final_interactions_Protein.List_Interaction_Protein{protein_table1});
    end
    end
end

%Copy List_of_final_interactions_Protein_60pc to use for other experiuments
List_of_final_interactions_Protein_60pc=List_of_final_interactions_Protein;

%Sort structure by Center/Gaussian write out table
[location2,location2]=sort([List_of_final_interactions_Protein.Center]); %Sort strcuture based on protein name
for sort_counter1=1:Maximum_observed_interactions(1)
Sorted_of_final_interactions_Center_60pc.Guassian{sort_counter1,1} = List_of_final_interactions_Protein.Guassian{location2(sort_counter1)};
Sorted_of_final_interactions_Center_60pc.Protein{sort_counter1,1}= List_of_final_interactions_Protein.Protein{location2(sort_counter1)};
Sorted_of_final_interactions_Center_60pc.Protein_NoIsoform{sort_counter1,1}=List_of_final_interactions_Protein.Protein_NoIsoform{location2(sort_counter1)};
    for MajorID_protein_counter = 1:dimension_Protein_MajorID_NoIsoforms(2)
    Sorted_of_final_interactions_Center_60pc.Protein_MajorID_NoIsoforms{sort_counter1,MajorID_protein_counter}=List_of_final_interactions_Protein.Protein_MajorID_NoIsoforms{location2(sort_counter1),MajorID_protein_counter};    
    end
Sorted_of_final_interactions_Center_60pc.Center{sort_counter1,1}= List_of_final_interactions_Protein.Center(location2(sort_counter1));
    for sort_counter2=1:Maximum_observed_interactions(2)
    Sorted_of_final_interactions_Center_60pc.Interaction_Gaussian{sort_counter1,sort_counter2}= List_of_final_interactions_Protein.Interaction_Gaussian{location2(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Center_60pc.Interaction_Protein{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Protein{location2(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Center_60pc.Interaction_Center{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Center(location2(sort_counter1),sort_counter2);
    Sorted_of_final_interactions_Center_60pc.Interaction_known_in_corum{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_known_in_corum(location2(sort_counter1),sort_counter2);
    end
Sorted_of_final_interactions_Center_60pc.Number_Interaction_known_in_corum(sort_counter1,1)= List_of_final_interactions_Protein.Number_Interaction_known_in_corum(location2(sort_counter1));
Sorted_of_final_interactions_Center_60pc.Total_Interaction(sort_counter1,1)=List_of_final_interactions_Protein.Total_Interaction(location2(sort_counter1));
Sorted_of_final_interactions_Center_60pc.List_Interaction_Gaussians(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Gaussians(location2(sort_counter1));
Sorted_of_final_interactions_Center_60pc.List_Interaction_Protein(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Protein(location2(sort_counter1));
end

%write out sorted structure array, Gaussian table
Sorted_Output_Interactions_60pc_precision_name=strcat('Sorted_Output_Interactions_Gaussians_60pc_replicate',mat2str(replicate_counter),'.csv');
fid_sorted_interactions_60_1 = fopen(Sorted_Output_Interactions_60pc_precision_name,'w');
fprintf (fid_sorted_interactions_60_1,'%s,%s,%s,%s,%s,%s,%s,\n',...
'Guassian', 'Protein_name(Isoform)','Protein_name(NoIsoform)','Center','Number of interactions','Number of interactions in Corum','Protein Interactions determined');               %Write Header
for write_out_sorted=1:Maximum_observed_interactions(1)
       fprintf(fid_sorted_interactions_60_1, '%s,%s,%s,%6.4f,%6.4f,%6.4f,%s,\n',...
       Sorted_of_final_interactions_Center_60pc.Guassian{write_out_sorted},...
       Sorted_of_final_interactions_Center_60pc.Protein{write_out_sorted},...
       Sorted_of_final_interactions_Center_60pc.Protein_NoIsoform{write_out_sorted},...
       Sorted_of_final_interactions_Center_60pc.Center{write_out_sorted},...
       Sorted_of_final_interactions_Center_60pc.Total_Interaction{write_out_sorted},...
       Sorted_of_final_interactions_Center_60pc.Number_Interaction_known_in_corum(write_out_sorted),...
       Sorted_of_final_interactions_Center_60pc.List_Interaction_Gaussians{write_out_sorted});
end
fclose(fid_sorted_interactions_60_1);

%Sort strcuture by Protein name
[location,location]=sort([List_of_final_interactions_Protein.Protein]); %Sort strcuture based on protein name
for sort_counter1=1:Maximum_observed_interactions(1)
Sorted_of_final_interactions_Protein_60pc.Guassian{sort_counter1,1} = List_of_final_interactions_Protein.Guassian{location(sort_counter1)};
Sorted_of_final_interactions_Protein_60pc.Protein{sort_counter1,1}= List_of_final_interactions_Protein.Protein{location(sort_counter1)};
Sorted_of_final_interactions_Protein_60pc.Protein_NoIsoform{sort_counter1,1}=List_of_final_interactions_Protein.Protein_NoIsoform{location(sort_counter1)};
    for MajorID_protein_counter = 1:dimension_Protein_MajorID_NoIsoforms(2)
    Sorted_of_final_interactions_Protein_60pc.Protein_MajorID_NoIsoforms{sort_counter1,MajorID_protein_counter}=List_of_final_interactions_Protein.Protein_MajorID_NoIsoforms{location(sort_counter1),MajorID_protein_counter};    
    end
Sorted_of_final_interactions_Protein_60pc.Center{sort_counter1,1}= List_of_final_interactions_Protein.Center(location(sort_counter1));
    for sort_counter2=1:Maximum_observed_interactions(2)
    Sorted_of_final_interactions_Protein_60pc.Interaction_Gaussian{sort_counter1,sort_counter2}= List_of_final_interactions_Protein.Interaction_Gaussian{location(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Protein_60pc.Interaction_Protein{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Protein{location(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Protein_60pc.Interaction_Center{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Center(location(sort_counter1),sort_counter2);
    Sorted_of_final_interactions_Protein_60pc.Interaction_known_in_corum{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_known_in_corum(location(sort_counter1),sort_counter2);
    end
Sorted_of_final_interactions_Protein_60pc.Number_Interaction_known_in_corum(sort_counter1,1)= List_of_final_interactions_Protein.Number_Interaction_known_in_corum(location(sort_counter1));
Sorted_of_final_interactions_Protein_60pc.Total_Interaction(sort_counter1,1)=List_of_final_interactions_Protein.Total_Interaction(location(sort_counter1));
Sorted_of_final_interactions_Protein_60pc.List_Interaction_Gaussians(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Gaussians(location(sort_counter1));
Sorted_of_final_interactions_Protein_60pc.List_Interaction_Protein(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Protein(location(sort_counter1));    
end

%write out sorted structure array, protein table
Sorted_Output_Interactions_Protein_60pc_name=strcat('Sorted_Output_Interactions_Proteins_60pc_replicate',mat2str(replicate_counter),'.csv');
fid_sorted_interactions_60_2 = fopen(Sorted_Output_Interactions_Protein_60pc_name,'w');
fprintf (fid_sorted_interactions_60_2,'%s,%s,%s,%s,%s,%s,%s,\n',...
'Guassian', 'Protein_name (Isoform)','Protein_name (NoIsoform)','Center','Number of interactions','Number of interactions in Corum','Protein Interactions determined');               %Write Header
for write_out_sorted=1:Maximum_observed_interactions(1)
       fprintf(fid_sorted_interactions_60_2, '%s,%s,%s,%6.4f,%6.4f,%6.4f,%s,\n',...
       Sorted_of_final_interactions_Protein_60pc.Guassian{write_out_sorted},...
       Sorted_of_final_interactions_Protein_60pc.Protein{write_out_sorted},...
       Sorted_of_final_interactions_Protein_60pc.Protein_NoIsoform{write_out_sorted},...
       Sorted_of_final_interactions_Protein_60pc.Center{write_out_sorted},...
       Sorted_of_final_interactions_Protein_60pc.Total_Interaction{write_out_sorted},...
       Sorted_of_final_interactions_Protein_60pc.Number_Interaction_known_in_corum(write_out_sorted),...
       Sorted_of_final_interactions_Protein_60pc.List_Interaction_Protein{write_out_sorted});
end
fclose(fid_sorted_interactions_60_2);

%% Define precision required 50pc
[~, idx] = min(abs(Precision_replicate_level - 0.5)); %determines the postion which gives a precision closest to 0.5
distance_with_precision_abover_0_5= test_precision(1,idx); %determines the point precision which gives a precision closest to 0.5

%note This module has been to designed to use the settings which give the highest number of high precision interactions with the widest Euc
Required_precision = distance_with_precision_abover_0_5; %set precision
Centre_Euc_with_require_precision = optimisation_matrix_precision > Required_precision;
Result_matrix =optimisation_matrix_number_interaction;
precision_counter= Optimisation_Dimension*Optimisation_Dimension; %determined by the number of parameters checked during the optimisation, thus if 200*200 matrix this number is 40000
parfor v=1:precision_counter   % merge the results of checking the precision and number of interactions
        if Centre_Euc_with_require_precision(v)==1
                Result_matrix(v) = optimisation_matrix_number_interaction(v);
                    else
                       Result_matrix(v)=NaN;
        end             
end
Average_interactions_columns=nanmean(Result_matrix); %Use to measure what the PPI landscape as Center varies with desired precision 
Average_interactions=nanmean(Average_interactions_columns); %Use to measure what the global PPI landscape over Euc and Center

optimisation_matrix_required_precision_name=strcat('optimisation_matrix_interaction_numbers_with_required_precision_50pc_replicate',mat2str(replicate_counter),'.csv');
csvwrite(optimisation_matrix_required_precision_name,Result_matrix); %write out the result matrix for use later

%determine the Euc and Center which give the user defined precision

maximum_counter=Optimisation_Dimension; %define based on the optimisation_matrix number 
Interactions_based_one_summed_precisions = zeros(Dimensions_Gaus_import,Dimensions_Gaus_import);
parfor vi=1:maximum_counter   % scan down the Euc row and find which values results the required precision and maximum number of interactions
    Maximum_interactions_at_given_Center = max(Result_matrix(vi,:));
    if ~isnan(Maximum_interactions_at_given_Center)
    inds =find(Result_matrix(vi,:)==Maximum_interactions_at_given_Center);
    [col] = inds(1); %this takes the first number of inds thus ensure the most conservative Center will be used
    optimised_center = Gaussian_fit_Distance_optimisation(vi);
    optimised_Euc = Euclidian_Distance_optimisation(col);
    
    %use optimisation code to maximum out
    MatrixGaussian = Dist.Gaussian_fits < optimised_center;
    Matrix_Distance = Dist.Euc < optimised_Euc;

    % determine final interaction matrices
    Interactions__precisions = (MatrixGaussian & MatrixHeight & MatrixWidth & Matrix_Distance & inverse_self);
    Interactions_based_one_summed_precisions = (Interactions_based_one_summed_precisions | Interactions__precisions);
    end
end

% determine final interaction matrices
Final_interactions = (Final_Interactions_Based_only_Eucldian_distance | Interactions_based_one_summed_precisions); %test just interactions on Euc dist

        % Calculate True possitives
        TP_Matrix_Int = (Final_interactions & TP_Matrix);
        TP= length(find(TP_Matrix_Int==1));

        % Calculate false negatives
        MaxTP= length(find(TP_Matrix==1));
        FN= MaxTP-TP;

        % Calculate false positives
        Int2_matrix= (Final_interactions & Int_matrix);
        Int3_matrix= (Int2_matrix & inverse_self);
        FP= length(find(Int3_matrix==1))-TP;
        
        % Find True negatives
        TN= (length(find(Int_matrix==1)))-FP;
        
        %Generate negative matrix for final interaction analysis
        Negative1_matrix=(Int_matrix  & ~Final_interactions);
        Negative2_matrix=(Negative1_matrix & inverse_self);        

%% final results
Recall = TP/(TP+FN)
Precision = TP/(TP+FP)
TPR = TP/(TP+FN)
FPR =FP/(FP+TN)
No_Int= length(find(triu(Final_interactions==1)))

Summary_Results_50pc_name=strcat('Summary_Results_50pc_replicate',mat2str(replicate_counter),'.csv');
fid_50_2= fopen(Summary_Results_50pc_name,'wt'); % create the summary file of the interaction output 
fprintf (fid_50_2,'%s,%s,%s,%s,%s,\n',... %header for OutputGaus output
                'Recall', 'Precision', 'TPR', 'FPR', 'No_Int'); %Write Header
fprintf (fid_50_2,'%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,\n',...
                    Recall, Precision, TPR, FPR, No_Int);
fprintf (fid_50_2,',\n');             
fclose(fid_50_2);

%% write out final interaction matrix for use to determine the global FDR, TPR, recall and precision across the experiments
Determined_interactions_50pc_precision_name=strcat('Determined_interactions_50pc_precision_replicate',mat2str(replicate_counter),'.csv');
csvwrite(Determined_interactions_50pc_precision_name,Final_interactions); %write out the result matrix for use later

%% Find "Binary" interactions

Protein_elution= strcat(num2str(C),'*',Protein.Isoform);
[n,m]=size(Dist.Euc);
Interactions =cell(1,1);

Protein_interactions=cell(1,2);
Protein_interactions_center=zeros(1,2);
Delta_EucDist=zeros(1,1);
Delta_Center=zeros(1,1);
Delta_Height=zeros(1,1);
Delta_Width=zeros(1,1);
Both_proteins_Corum=zeros(1,1);
Interaction_in_Corum=zeros(1,1);
final_count=0;
U = triu(Final_interactions);    %Take half of matrix therefor unique interactions
for vii=1:n
    for viii= 1:m
        if U(vii,viii)==1
            Int1=Protein_elution(vii);
            Int2=Protein_elution(viii);
            if strcmp(Int1, Int2)
            else
                final_count=1+final_count;
                Interactions{final_count,1}= strcat(Int1,'-',Int2);
                Delta_EucDist(final_count,1)= Dist.Euc(vii,viii);
                Delta_Center(final_count,1)= abs(C(vii)-C(viii));
                Delta_Height(final_count,1)= abs(H(vii)-H(viii));
                Delta_Width(final_count,1)= abs(W(vii)-W(viii));
                Both_proteins_Corum(final_count,1)= Int_matrix(vii,viii);
                Interaction_in_Corum(final_count,1)= TP_Matrix(vii,viii);
                Protein_interactions(final_count,1)=Protein.Isoform(vii);
                Protein_interactions(final_count,2)=Protein.Isoform(viii);
                Protein_interactions_center(final_count,1)=C(vii);
                Protein_interactions_center(final_count,2)=C(viii);
            end
        else
        end   
        
    end   
end

Num_Inter=length(Interactions);
QQ=([Interactions{:}])';

%Print interactions
Output_Interactions_50pc_precision_name=strcat('Output_Interactions_50pc_precision_replicate',mat2str(replicate_counter),'.csv');
fid11 = fopen(Output_Interactions_50pc_precision_name,'w');
fprintf (fid11,'%s,%s,%s,%s,%s,%s,%s,\n',...
        'Interaction', 'Delta Height', 'Delta Center', 'Delta Width', 'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for ix=1:Num_Inter    
fprintf(fid11,'%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',QQ{ix}, Delta_Height(ix), Delta_Center(ix), Delta_Width(ix), Delta_EucDist(ix), Both_proteins_Corum(ix), Interaction_in_Corum(ix));
end
fclose(fid11);

List_Interactions_50pc_precision_name=strcat('Interactions_list_50pc_precision_replicate',mat2str(replicate_counter),'.csv');
fid11B = fopen(List_Interactions_50pc_precision_name,'w');
fprintf (fid11B,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
        'Replicate','Protein A','Protein B', 'Center A', 'Center B',...
        'Delta Height', 'Delta Center', 'Delta Width',...
        'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for ix=1:Num_Inter    
fprintf(fid11B,'%6.3f,%s,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',replicate_counter,...
    Protein_interactions{ix,1},Protein_interactions{ix,2},Protein_interactions_center(ix,1),Protein_interactions_center(ix,2),...
    Delta_Height(ix),Delta_Center(ix), Delta_Width(ix), Delta_EucDist(ix),...
    Both_proteins_Corum(ix), Interaction_in_Corum(ix));
end
fclose(fid11B);

%% Create list of predicted negative interactions (FN +TN)
%Determine the false "binary" interaction list
U1=triu(Negative2_matrix); %Take half of matrix therefor unique interactions

%change directory
cd('../');
name_negative_file=strcat('Negative interactions_50pc_replicate_',mat2str(replicate_counter));
mkdir(name_negative_file);
parfor vii=1:n
    for viii= 1:m
        if U1(vii,viii)==1
            Int1=Protein_elution(vii);
            Int2=Protein_elution(viii);
            if strcmp(Int1, Int2)
            else
                FN_Interactions= {strcat(Int1,'-',Int2)};
                FN_Delta_EucDist= Dist.Euc(vii,viii);
                FN_Delta_Center= abs(C(vii)-C(viii));
                FN_Delta_Height= abs(H(vii)-H(viii));
                FN_Delta_Width= abs(W(vii)-W(viii));
                FN_Both_proteins_Corum= Int_matrix(vii,viii);
                FN_Interaction_in_Corum= TP_Matrix(vii,viii);
                FN_Protein_interactions1=Protein.Isoform(vii);
                FN_Protein_interactions2=Protein.Isoform(viii);
                FN_Protein_interactions_center1=C(vii);
                FN_Protein_interactions_center2=C(viii);
                
                %Write out data file
                fileName = strcat(mat2str(vii),'_',mat2str(viii),'_Interaction','.csv');
                fid1A = fopen(fileName,'at');
                fprintf(fid1A,'%6f,%s,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',...
                replicate_counter,FN_Protein_interactions1{1},...
                FN_Protein_interactions2{1},FN_Protein_interactions_center1(1),...
                FN_Protein_interactions_center2(1), FN_Delta_Height(1), FN_Delta_Center(1),...
                FN_Delta_Width(1),FN_Delta_EucDist(1), FN_Both_proteins_Corum(1), FN_Interaction_in_Corum(1));
                fclose(fid1A);
                movefile(fileName, name_negative_file);
            end
        else
        end   
        
    end   
end


%read dictory and concatinate negative interactions into a single file
dirData2 = dir(name_negative_file);      %# Get directory contents
dirData2 = dirData2(~[dirData2.isdir]);  %# Use only the file data
negative_fileNames = {dirData2.name};           %# Get file names

ExpName1 = strcat('^(\d+)\_(\d+)\','_Interaction','.csv');
Incorum_Neg_fileData = regexp(negative_fileNames,ExpName1,'tokens');  %# Find tokens
index1 = ~cellfun('isempty',Incorum_Neg_fileData);        %# Find index of matches
negative_fileNames = negative_fileNames(index1);        %# remove empty cells
Incorum_Neg_fileData = [Incorum_Neg_fileData{index1}];                %# Remove non-matches
Incorum_Neg_fileData = vertcat(Incorum_Neg_fileData{:});             %# Format token data
              
%sorting protein interactions in Corum in correct numerical order
nFiles1 = size(Incorum_Neg_fileData,1);              %# Number of files matching format
Incorum_Neg_fileData = str2double(Incorum_Neg_fileData);        %# Convert from strings to numbers
[Incorum_Neg_fileData,index1] = sortrows(Incorum_Neg_fileData,1:2);  %# Sort numeric values
negative_fileNames =  negative_fileNames(index1);               %# Apply sort to file names

Number_of_Neg_interactions = length(negative_fileNames);
Combined_Neg_interactions_data=cell(Number_of_Neg_interactions,11);

for xiv = 1:Number_of_Neg_interactions 
    %load file of interest
    fileName7 = strcat(name_negative_file,'\',negative_fileNames(xiv));
    File_to_be_read = fileName7{1};
    
    %import data files
    f = fopen (File_to_be_read);
    Processed_data = textscan(f, '%s %s %s %s %s %s %s %s %s %s %s %s %s', 'Delimiter', ',');
    fclose(f);  

    %Copy values to arrays to use in downstream analysis
    %Value 1, As the first value is the replicate there is not need to save
    %value
    
    %Value 2
    prot1_protein_names_entry=Processed_data{2};
    Neg_Corum_prot1_protein_names{xiv}=prot1_protein_names_entry{1};
    %Value 3
    prot2_protein_names_entry=Processed_data{3};
    Neg_Corum_prot2_protein_names{xiv}=prot2_protein_names_entry{1};
    %Value 4
    prot1_center=Processed_data{4};
    Neg_Corum_prot1_center(xiv)=str2double(prot1_center(1));
    %Value 5
    prot2_center=Processed_data{5};
    Neg_Corum_prot2_center(xiv)=str2double(prot2_center{1});
    %Value 6
    Neg_Delta_height_value=Processed_data{6};
    Neg_Delta_height(xiv)=str2double(Neg_Delta_height_value{1});
    %Value 7
    Neg_Center_value=Processed_data{7};
    Neg_Center(xiv)=str2double(Neg_Center_value(1));
    %Value 8
    Neg_Width_value=Processed_data{8};
    Neg_Width(xiv)=str2double(Neg_Width_value(1));
    %Value 9
    Neg_Euc_value=Processed_data{9};
    Neg_Euc(xiv)=str2double(Neg_Euc_value{1});
    %Value 10
    Neg_both_in_corum_value=Processed_data{10};
    Neg_both_in_corum(xiv)=str2double(Neg_both_in_corum_value{1});
    %Value 11
    Neg_interaction_in_corum_value=Processed_data{11};
    Neg_interaction_in_corum(xiv)=str2double(Neg_interaction_in_corum_value(1));
end

%Write out negative interaction list, used to create the final TPR, FPR
%precision and recall values
cd('ROC analysis data');
Neg_proteins_interactions_name=strcat('Neg_proteins_interactions_50pc_replicate_',mat2str(replicate_counter),'.csv');
fid_Neg= fopen(Neg_proteins_interactions_name,'w'); 
fprintf (fid_Neg,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
        'Replicate','Protein A','Protein B', 'Center A', 'Center B',...
        'Delta Height', 'Delta Center', 'Delta Width',...
        'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for xv=1:Number_of_Neg_interactions   
fprintf(fid_Neg,'%6f,%s,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',...
replicate_counter(1) ,Neg_Corum_prot1_protein_names{xv}, Neg_Corum_prot2_protein_names{xv},...
Neg_Corum_prot1_center(xv),Neg_Corum_prot2_center(xv), Neg_Delta_height(xv),...
Neg_Center(xv),Neg_Width(xv),Neg_Euc(xv),Neg_both_in_corum(xv),Neg_interaction_in_corum(xv));
end     
fclose(fid_Neg);

%% List of interactions
List_of_final_interactions = Final_interactions;    %As we want the complete list do not use triu
Protein_elution_1= strcat(num2str(C),'*',Protein.Isoform);
[n,m]=size(Dist.Euc);
row_counter=0;
empty_cell=cell(1,1);
for vii=1:n
    List_of_final_interactions_Protein.Guassian(vii,1)= Protein_elution_1(vii);
    List_of_final_interactions_Protein.Protein(vii,1)= Protein.Isoform(vii);
    List_of_final_interactions_Protein.Protein_NoIsoform(vii,1)= Protein.NoIsoform(vii);
        for MajorID_protein_counter = 1:dimension_Protein_MajorID_NoIsoforms(2) 
        List_of_final_interactions_Protein.Protein_MajorID_NoIsoforms(vii,MajorID_protein_counter)= Protein.MajorID_NoIsoforms(vii,MajorID_protein_counter);
        end
    List_of_final_interactions_Protein.Center(vii,1)= C(vii);
    final_count=0;
    for viii= 1:m
        if List_of_final_interactions(vii,viii)== 1 %Determine if proteins interact
            final_count=1+final_count;
            List_of_final_interactions_Protein.Interaction_Gaussian(vii,final_count)=Protein_elution_1(viii);
            List_of_final_interactions_Protein.Interaction_Protein(vii,final_count)=Protein.Isoform(viii);
            List_of_final_interactions_Protein.Interaction_Protein_NoIsoform(vii,final_count)=Protein.NoIsoform(viii);
            List_of_final_interactions_Protein.Interaction_Center(vii,final_count)=C(viii);
            if  TP_Matrix(vii,viii)== 1 %Determine if interaction is known in Corum
                List_of_final_interactions_Protein.Interaction_known_in_corum(vii,final_count)= 1;
            else
                List_of_final_interactions_Protein.Interaction_known_in_corum(vii,final_count)= 0;
            end
        end
    end
    if final_count == 0
    List_of_final_interactions_Protein.Interaction_Gaussian(vii,1)= {NaN};
    List_of_final_interactions_Protein.Interaction_Protein(vii,1)={NaN};
    List_of_final_interactions_Protein.Interaction_Center(vii,1)= 0;
    List_of_final_interactions_Protein.Interaction_known_in_corum(vii,1)=0;
    end
end

%Remove values in List of final with NaN
List_of_final_interactions_Protein.Interaction_Gaussian(cellfun(@(x) any(isnan(x)),List_of_final_interactions_Protein.Interaction_Gaussian)) = {''};
List_of_final_interactions_Protein.Interaction_Protein(cellfun(@(x) any(isnan(x)),List_of_final_interactions_Protein.Interaction_Protein)) = {''};

%Determine the dimension of List_of_final_interactions_Protein.Interaction_Protein array
Maximum_observed_interactions=size(List_of_final_interactions_Protein.Interaction_Protein); % determine the dimension of the interaction array in structure

%Count number of corum interactions and total interactions
List_of_final_interactions_Protein.List_Interaction_Gaussians=cell(Maximum_observed_interactions(1),1);
List_of_final_interactions_Protein.Total_Interaction=cell(Maximum_observed_interactions(1),1);
for Interaction_number_counter1=1:Maximum_observed_interactions(1)
    Interaction_counter_for_protein=0;
    for Interaction_number_counter2=1:Maximum_observed_interactions(2)
    Guassian_to_combine(Interaction_number_counter2)=List_of_final_interactions_Protein.Interaction_Gaussian(Interaction_number_counter1,Interaction_number_counter2);
    if any(Guassian_to_combine{Interaction_number_counter2})== 1 % test if value is not zero
        String_name= strcat(Guassian_to_combine{Interaction_number_counter2},' ; ');
        List_of_final_interactions_Protein.List_Interaction_Gaussians{Interaction_number_counter1}=strcat(String_name,List_of_final_interactions_Protein.List_Interaction_Gaussians{Interaction_number_counter1});
        Interaction_counter_for_protein=Interaction_counter_for_protein+1;
    end
    end
    List_of_final_interactions_Protein.Total_Interaction{Interaction_number_counter1,1}=Interaction_counter_for_protein;
end

%Sum number of known interactions in corum for protein
for Interaction_number_counter1=1:Maximum_observed_interactions(1)
for Interaction_number_counter2=1:Maximum_observed_interactions(2)
Detected_Center_Interaction_known_in_corum(Interaction_number_counter2)= List_of_final_interactions_Protein.Interaction_known_in_corum(Interaction_number_counter1,Interaction_number_counter2);
end
List_of_final_interactions_Protein.Number_Interaction_known_in_corum(Interaction_number_counter1,1)= sum(Detected_Center_Interaction_known_in_corum);
end

%Protein list of interactions
List_of_final_interactions_Protein.List_Interaction_Protein=cell(Maximum_observed_interactions(1),1);
for protein_table1=1:Maximum_observed_interactions(1)
    for protein_table2=1:Maximum_observed_interactions(2)
    Proteins_to_combine(protein_table2)=List_of_final_interactions_Protein.Interaction_Protein(protein_table1,protein_table2);
    if any(Proteins_to_combine{protein_table2})== 1 % test if value is not zero
        String_name= strcat(Proteins_to_combine{protein_table2},' ; ');
        List_of_final_interactions_Protein.List_Interaction_Protein{protein_table1,1}=strcat(String_name,List_of_final_interactions_Protein.List_Interaction_Protein{protein_table1});
    end
    end
end

%Copy List_of_final_interactions_Protein_50pc to use for other experiuments
List_of_final_interactions_Protein_50pc=List_of_final_interactions_Protein;

%Sort structure by Center/Gaussian write out table
[location2,location2]=sort([List_of_final_interactions_Protein.Center]); %Sort strcuture based on protein name
for sort_counter1=1:Maximum_observed_interactions(1)
Sorted_of_final_interactions_Center_50pc.Guassian{sort_counter1,1} = List_of_final_interactions_Protein.Guassian{location2(sort_counter1)};
Sorted_of_final_interactions_Center_50pc.Protein{sort_counter1,1}= List_of_final_interactions_Protein.Protein{location2(sort_counter1)};
Sorted_of_final_interactions_Center_50pc.Protein_NoIsoform{sort_counter1,1}=List_of_final_interactions_Protein.Protein_NoIsoform{location2(sort_counter1)};
    for MajorID_protein_counter = 1:dimension_Protein_MajorID_NoIsoforms(2)
    Sorted_of_final_interactions_Center_50pc.Protein_MajorID_NoIsoforms{sort_counter1,MajorID_protein_counter}=List_of_final_interactions_Protein.Protein_MajorID_NoIsoforms{location2(sort_counter1),MajorID_protein_counter};    
    end
Sorted_of_final_interactions_Center_50pc.Center{sort_counter1,1}= List_of_final_interactions_Protein.Center(location2(sort_counter1));
    for sort_counter2=1:Maximum_observed_interactions(2)
    Sorted_of_final_interactions_Center_50pc.Interaction_Gaussian{sort_counter1,sort_counter2}= List_of_final_interactions_Protein.Interaction_Gaussian{location2(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Center_50pc.Interaction_Protein{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Protein{location2(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Center_50pc.Interaction_Center{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Center(location2(sort_counter1),sort_counter2);
    Sorted_of_final_interactions_Center_50pc.Interaction_known_in_corum{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_known_in_corum(location2(sort_counter1),sort_counter2);
    end
Sorted_of_final_interactions_Center_50pc.Number_Interaction_known_in_corum(sort_counter1,1)= List_of_final_interactions_Protein.Number_Interaction_known_in_corum(location2(sort_counter1));
Sorted_of_final_interactions_Center_50pc.Total_Interaction(sort_counter1,1)=List_of_final_interactions_Protein.Total_Interaction(location2(sort_counter1));
Sorted_of_final_interactions_Center_50pc.List_Interaction_Gaussians(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Gaussians(location2(sort_counter1));
Sorted_of_final_interactions_Center_50pc.List_Interaction_Protein(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Protein(location2(sort_counter1));
end

%write out sorted structure array, Gaussian table
Sorted_output_Interactions_Gaussian_50pc_name=strcat('Sorted_Output_Interactions_Gaussians_50pc_replicate',mat2str(replicate_counter),'.csv');
fid_sorted_interactions_50_1 = fopen(Sorted_output_Interactions_Gaussian_50pc_name,'w');
fprintf (fid_sorted_interactions_50_1,'%s,%s,%s,%s,%s,%s,%s,\n',...
'Guassian', 'Protein_name(Isoform)','Protein_name(NoIsoform)','Center','Number of interactions','Number of interactions in Corum','Protein Interactions determined');               %Write Header
for write_out_sorted=1:Maximum_observed_interactions(1)
       fprintf(fid_sorted_interactions_50_1, '%s,%s,%s,%6.4f,%6.4f,%6.4f,%s,\n',...
       Sorted_of_final_interactions_Center_50pc.Guassian{write_out_sorted},...
       Sorted_of_final_interactions_Center_50pc.Protein{write_out_sorted},...
       Sorted_of_final_interactions_Center_50pc.Protein_NoIsoform{write_out_sorted},...
       Sorted_of_final_interactions_Center_50pc.Center{write_out_sorted},...
       Sorted_of_final_interactions_Center_50pc.Total_Interaction{write_out_sorted},...
       Sorted_of_final_interactions_Center_50pc.Number_Interaction_known_in_corum(write_out_sorted),...
       Sorted_of_final_interactions_Center_50pc.List_Interaction_Gaussians{write_out_sorted});
end
fclose(fid_sorted_interactions_50_1);

%Sort strcuture by Protein name
[location,location]=sort([List_of_final_interactions_Protein.Protein]); %Sort strcuture based on protein name
for sort_counter1=1:Maximum_observed_interactions(1)
Sorted_of_final_interactions_Protein_50pc.Guassian{sort_counter1,1} = List_of_final_interactions_Protein.Guassian{location(sort_counter1)};
Sorted_of_final_interactions_Protein_50pc.Protein{sort_counter1,1}= List_of_final_interactions_Protein.Protein{location(sort_counter1)};
Sorted_of_final_interactions_Protein_50pc.Protein_NoIsoform{sort_counter1,1}=List_of_final_interactions_Protein.Protein_NoIsoform{location(sort_counter1)};
    for MajorID_protein_counter = 1:dimension_Protein_MajorID_NoIsoforms(2)
    Sorted_of_final_interactions_Protein_50pc.Protein_MajorID_NoIsoforms{sort_counter1,MajorID_protein_counter}=List_of_final_interactions_Protein.Protein_MajorID_NoIsoforms{location(sort_counter1),MajorID_protein_counter};    
    end
Sorted_of_final_interactions_Protein_50pc.Center{sort_counter1,1}= List_of_final_interactions_Protein.Center(location(sort_counter1));
    for sort_counter2=1:Maximum_observed_interactions(2)
    Sorted_of_final_interactions_Protein_50pc.Interaction_Gaussian{sort_counter1,sort_counter2}= List_of_final_interactions_Protein.Interaction_Gaussian{location(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Protein_50pc.Interaction_Protein{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Protein{location(sort_counter1),sort_counter2};
    Sorted_of_final_interactions_Protein_50pc.Interaction_Center{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_Center(location(sort_counter1),sort_counter2);
    Sorted_of_final_interactions_Protein_50pc.Interaction_known_in_corum{sort_counter1,sort_counter2}=List_of_final_interactions_Protein.Interaction_known_in_corum(location(sort_counter1),sort_counter2);
    end
Sorted_of_final_interactions_Protein_50pc.Number_Interaction_known_in_corum(sort_counter1,1)= List_of_final_interactions_Protein.Number_Interaction_known_in_corum(location(sort_counter1));
Sorted_of_final_interactions_Protein_50pc.Total_Interaction(sort_counter1,1)=List_of_final_interactions_Protein.Total_Interaction(location(sort_counter1));
Sorted_of_final_interactions_Protein_50pc.List_Interaction_Gaussians(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Gaussians(location(sort_counter1));
Sorted_of_final_interactions_Protein_50pc.List_Interaction_Protein(sort_counter1,1)=List_of_final_interactions_Protein.List_Interaction_Protein(location(sort_counter1));    
end

%write out sorted structure array, protein table
Sorted_output_Interactions_Protein_50pc_name=strcat('Sorted_Output_Interactions_Protein_50pc_replicate',mat2str(replicate_counter),'.csv');
fid_sorted_interactions_50_2 = fopen(Sorted_output_Interactions_Protein_50pc_name,'w');
fprintf (fid_sorted_interactions_50_2,'%s,%s,%s,%s,%s,%s,%s,\n',...
'Guassian', 'Protein_name (Isoform)','Protein_name (NoIsoform)','Center','Number of interactions','Number of interactions in Corum','Protein Interactions determined');               %Write Header
for write_out_sorted=1:Maximum_observed_interactions(1)
       fprintf(fid_sorted_interactions_50_2, '%s,%s,%s,%6.4f,%6.4f,%6.4f,%s,\n',...
       Sorted_of_final_interactions_Protein_50pc.Guassian{write_out_sorted},...
       Sorted_of_final_interactions_Protein_50pc.Protein{write_out_sorted},...
       Sorted_of_final_interactions_Protein_50pc.Protein_NoIsoform{write_out_sorted},...
       Sorted_of_final_interactions_Protein_50pc.Center{write_out_sorted},...
       Sorted_of_final_interactions_Protein_50pc.Total_Interaction{write_out_sorted},...
       Sorted_of_final_interactions_Protein_50pc.Number_Interaction_known_in_corum(write_out_sorted),...
       Sorted_of_final_interactions_Protein_50pc.List_Interaction_Protein{write_out_sorted});
end
fclose(fid_sorted_interactions_50_2);

cd('../');

%%Comparsion of properties of Corum compared to network observed using
%%PCP-SEC-SILAC
cd('ROC analysis related Figures');

%Count the number of interactions found for each protein in corum
number_corum_interactions=zeros(number_of_Unique_Corum,1);
for corum_counter1=1:number_of_Unique_Corum
    number_corum_interactions(corum_counter1)=sum(strcmp(Unique_Corum(corum_counter1), Corum_Protein_names(:,1)));
end
    
%Count the number of interactions found for each protein in corum in our
%dataset
number_corum_interactions_in_our_dataset=zeros(number_of_pos_corum_proteins,1);
for corum_counter1=1:number_of_pos_corum_proteins
    number_corum_interactions_in_our_dataset(corum_counter1)=sum(strcmp(Pos_Corum_proteins(corum_counter1), Corum_Protein_names(:,1)));
end    

%Count the number of interactions found for each protein in corum in our
%dataset, for individual proteins

%column one total interactions, column two interactions also in corum
number_interactions_in_our_dataset_of_corum_protien_70pc=zeros(number_of_pos_corum_proteins,2);
number_interactions_in_our_dataset_of_corum_protien_60pc=zeros(number_of_pos_corum_proteins,2);
number_interactions_in_our_dataset_of_corum_protien_50pc=zeros(number_of_pos_corum_proteins,2);
%Array to store Center, column 1 isoform name, column 2 noisoform name
Name_of_protein_to_write_out=cell(number_of_pos_corum_proteins,2);

for corum_counter1=1:number_of_pos_corum_proteins
    [location_of_corum_protien]=ind2sub(Dimensions_Gaus_import(1), strmatch(Pos_Corum_proteins(corum_counter1), [Protein_location_array.MajorID_NoIsoforms{:}], 'exact'));
    row_numbers1=Protein_location_array.row(location_of_corum_protien);
    row_numbers2=unique([row_numbers1{:}]);
    number_of_rows=length(row_numbers2);
    %Record center
    Name_of_protein_to_write_out{corum_counter1,1}=Protein.NoIsoform(row_numbers2(1));
    Name_of_protein_to_write_out{corum_counter1,2}=Protein.Isoform(row_numbers2(1));
    %Create counter to store total interaction values for a given protein
    NES_70pc =0;
    NES_60pc =0;
    NES_50pc =0;
    %Create counter to store interaction in corum values for a given protein
    NES1_70pc =0;
    NES1_60pc =0;
    NES1_50pc =0;    
    %Check is multiple gaussians have been detected for a given protein
        for write_out_counter=1:number_of_rows
        %for interaction with precision >70%
        NES_70pc= NES_70pc + List_of_final_interactions_Protein_70pc.Number_Interaction_known_in_corum(row_numbers2(write_out_counter));
        NES1_70pc= NES1_70pc +  List_of_final_interactions_Protein_70pc.Total_Interaction{row_numbers2(write_out_counter)};
        %for interaction with precision >60%
        NES_60pc= NES_60pc + List_of_final_interactions_Protein_60pc.Number_Interaction_known_in_corum(row_numbers2(write_out_counter));
        NES1_60pc= NES1_60pc + List_of_final_interactions_Protein_60pc.Total_Interaction{row_numbers2(write_out_counter)};
        %for interaction with precision >50%
        NES_50pc= NES_50pc + List_of_final_interactions_Protein_50pc.Number_Interaction_known_in_corum(row_numbers2(write_out_counter));
        NES1_50pc= NES1_50pc + List_of_final_interactions_Protein_50pc.Total_Interaction{row_numbers2(write_out_counter)};
        end
    number_interactions_in_our_dataset_of_corum_protien_70pc(corum_counter1,1)= NES_70pc;
    number_interactions_in_our_dataset_of_corum_protien_60pc(corum_counter1,1)= NES_60pc;
    number_interactions_in_our_dataset_of_corum_protien_50pc(corum_counter1,1)= NES_50pc;
    number_interactions_in_our_dataset_of_corum_protien_70pc(corum_counter1,2)= NES1_70pc;
    number_interactions_in_our_dataset_of_corum_protien_60pc(corum_counter1,2)= NES1_60pc;
    number_interactions_in_our_dataset_of_corum_protien_50pc(corum_counter1,2)= NES1_50pc;
end  

%average number of interactions
average_70pc(replicate_counter)=mean([List_of_final_interactions_Protein_70pc.Total_Interaction{:}]);
average_60pc(replicate_counter)=mean([List_of_final_interactions_Protein_60pc.Total_Interaction{:}]);
average_50pc(replicate_counter)=mean([List_of_final_interactions_Protein_50pc.Total_Interaction{:}]);
average_interaction_corum(replicate_counter)=mean(number_corum_interactions);

%Write out summary of average number of interactions

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


%plot Interaction data
f1=figure
%plot1
subplot(4,1,1);
xvalues1 = 0:70;
[n1, xout1] = hist(number_interactions_in_our_dataset_of_corum_protien_70pc(:,2),xvalues1);
P1=bar(xout1,n1);
set(P1, 'FaceColor', colour_to_use(1,:)) 
xlim([0,60]);
ylim([0,700]);
title('Interactions Observed with precision >70%');
xlabel('Number of interaction pre protein');
ylabel('# of proteins');
set(gca,'YScale','log');
hold on
[n2, xout2] = hist(number_interactions_in_our_dataset_of_corum_protien_70pc(:,1),xvalues1);
P2= bar(xout2,n2);
set(P2, 'FaceColor', colour_to_use(4,:)) 
legend([P1 P2],{'Total interactions', 'Interactions observed in Corum'});
%plot2
subplot(4,1,2);
[n1, xout1] = hist(number_interactions_in_our_dataset_of_corum_protien_60pc(:,2),xvalues1);
P1=bar(xout1,n1);
set(P1, 'FaceColor', colour_to_use(1,:)) 
xlim([0,60]);
ylim([0,700]);
title('Interactions Observed with precision >60%');
xlabel('Number of interaction pre protein');
ylabel('# of proteins');
hold on
[n2, xout2] = hist(number_interactions_in_our_dataset_of_corum_protien_60pc(:,1),xvalues1);
P2= bar(xout2,n2);
set(P2, 'FaceColor', colour_to_use(4,:)); 
set(gca,'YScale','log');
%plot3
subplot(4,1,3);
[n1, xout1] = hist(number_interactions_in_our_dataset_of_corum_protien_50pc(:,2),xvalues1);
P1=bar(xout1,n1);
set(P1, 'FaceColor', colour_to_use(1,:)); 
xlim([0,60]);
ylim([0,700]);
title('Interactions Observed with precision >50%');
xlabel('Number of interaction pre protein');
ylabel('# of proteins');
hold on
[n2, xout2] = hist(number_interactions_in_our_dataset_of_corum_protien_50pc(:,1),xvalues1);
P2= bar(xout2,n2);
set(P2, 'FaceColor', colour_to_use(4,:));
set(gca,'YScale','log');
%plot4
subplot(4,1,4);
[n1, xout1] = hist(number_corum_interactions(:,1),xvalues1);
P1=bar(xout1,n1);
set(P1, 'FaceColor', colour_to_use(10,:));
xlim([0,60]);
ylim([0,700]);
title('Interactions Observed in corum');
xlabel('Number of interaction pre protein');
ylabel('# of proteins');
set(gca,'YScale','log');
Interactions_pre_proteins_name=strcat('Interactions_pre_proteins_replicate',mat2str(replicate_counter),'.png');
saveas(f1, Interactions_pre_proteins_name);

cd('../');
cd('../');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%make directory
mkdir('Combined results');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Write summary of analysis of from replicates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cd('Combined results');
fid_Eu2= fopen('Summary_Interactions_determined_based_on_Euclidean_distance.csv','wt'); % Eu dis interaction output 
fprintf (fid_Eu2,'%s,%s,%s,\n',... %header for output
'Replicate','Total number of interactions', 'Unique number of interactions'); %Write Header 
for replicate_counter=1:(number_of_replicates*number_of_channels)
fprintf(fid_Eu2,'%6f,%6.3f,%6.3f\n',...
    replicate_counter, Total_Interactions_Based_only_Eucldian_distance(replicate_counter), Count_Interactions_Based_only_Eucldian_distance(replicate_counter)); 
end
fclose(fid_Eu2);
cd('../');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Combine data to generate final interation networks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Generate combined list of all binary interactions (including Center
% information)

%generate empty table to place binary interactions
binary_interaction_header ={'Replicate','Protein A','Protein B', 'Center A', 'Center B',...
        'Delta Height', 'Delta Center', 'Delta Width',...
        'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'}; %Write Header

%Create precision values to use    
Precision_values=[70;60;50];
length_precision_values=3;

for precision_write_out_counter=1:length_precision_values


%Generate list of binary interactions
binary_interaction_list=binary_interaction_header;

%Generate list of all FN binary interactions
FN_binary_interaction_list=binary_interaction_header;

for replicate_counter=1:(number_of_replicates*number_of_channels)
    %Go to replicate directory
    cd(Save_replicate_folder_name{replicate_counter});
    
    %Read in list outputs to create final interaction network
    import_name=strcat('ROC analysis data\Interactions_list_',mat2str(Precision_values(precision_write_out_counter)),'pc_precision_replicate',mat2str(replicate_counter),'.csv');
      
    %Read in list outputs to create final interaction network
    FN_import_name=strcat('ROC analysis data\Neg_proteins_interactions_',mat2str(Precision_values(precision_write_out_counter)),'pc_replicate_',mat2str(replicate_counter),'.csv');

    %Import data
    import = fopen (import_name);
    import_data = textscan(import, '%s', 'Delimiter',',');
    fclose(import);
    
    %Import Negativedata
    FN_import = fopen (FN_import_name);
    FN_import_data = textscan(FN_import, '%s', 'Delimiter',',');
    fclose(FN_import);
    
    %Reshape data for use in analysis
    import_Dimension=size(import_data{:});
    import_reshaped=reshape(import_data{:},11,(import_Dimension(1)/11))';
    import_reshaped_size=size(import_reshaped);
    %Remove title
    import_reshaped=import_reshaped(2:end,:);
    
    %Process Negative interactions
    import_FN_Dimension=size(FN_import_data{:});
    import_FN_reshaped=reshape(FN_import_data{:},11,(import_FN_Dimension(1)/11))';
    import_FN_reshaped_size=size(import_FN_reshaped);
    %Remove title
    import_FN_reshaped=import_FN_reshaped(2:end,:);

    %Concatenate arrays vertically
    %True interactions
    binary_interaction_list = vertcat(binary_interaction_list,import_reshaped);
    %False interactions
    FN_binary_interaction_list = vertcat(FN_binary_interaction_list,import_FN_reshaped);
    cd('../');
end

%% create a list of all the unique interactions (within two fractions) and how of end they are observed across replicates
%create unique identifier for true interactions
interaction_pairs=strcat(binary_interaction_list(:,2),'_',binary_interaction_list(:,3));
unique_interaction = unique(interaction_pairs(2:end));
length_unique_inter=length(unique_interaction);

%Create array to store values
interaction_final.unique_interactions=cell(0,0);
interaction_final.replicate_numbers=cell(0,0);
interaction_final.proteinA=cell(0,0);
interaction_final.proteinB=cell(0,0);
interaction_final.CenterA=cell(0,0);
interaction_final.CenterB=cell(0,0);
interaction_final.DeltaHeight=cell(0,0);
interaction_final.DeltaCenter=cell(0,0);
interaction_final.DeltaWidth=cell(0,0);
interaction_final.DeltaEucDist=cell(0,0);   
interaction_final.proteinInCorum=cell(0,0);
interaction_final.interactionInCorum=cell(0,0);

Unique_interaction_counter=1;
for interaction_counter2A=90:length_unique_inter
%Find location of unique protein interaction    
[location_interaction_pairs]=ind2sub(length_unique_inter, strmatch(unique_interaction(interaction_counter2A), [interaction_pairs(2:end)], 'exact'));

%find which replicate the interactions were seen in
while ~isempty(location_interaction_pairs(:))==1
    
   %determine if the centers of the interactions are withing 2 fractions
   num_interaction_pairs=length(location_interaction_pairs);
   center_interaction_pairs=binary_interaction_list(location_interaction_pairs+1,4);
   %test are center values within two
   Values_to_test=(abs(cellfun(@str2num,center_interaction_pairs)...
       -cellfun(@str2num,center_interaction_pairs(num_interaction_pairs(1)))))<2;
   %find positions within two
   position_within_two=find(Values_to_test == 1);
   
   %Define row counter
   row_counter=1;
   %Check if multiple centers are withing two fractions of each other
   if ~nnz(position_within_two)==0
   for counter=1:length(position_within_two)
   interaction_final.unique_interactions(Unique_interaction_counter,1)= unique_interaction(interaction_counter2A);
   %check size of replicate_numbers
   interaction_final.replicate_numbers{Unique_interaction_counter,row_counter}=...
   cellfun(@str2num,binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,1));
   % copy the values for the first unique interaction to the list
   interaction_final.proteinA(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,2);
   interaction_final.proteinB(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,3);
   interaction_final.CenterA(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,4);
   interaction_final.CenterB(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,5);
   interaction_final.DeltaHeight(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,6);
   interaction_final.DeltaCenter(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,7);
   interaction_final.DeltaWidth(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,8);
   interaction_final.DeltaEucDist(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,9);   
   interaction_final.proteinInCorum(Unique_interaction_counter,1)= binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,10);
   interaction_final.interactionInCorum(Unique_interaction_counter,1)= binary_interaction_list(location_interaction_pairs(position_within_two(counter))+1,11);
   
   %Add one to counters and remove location already checked
   row_counter=row_counter+1;
   end
   location_interaction_pairs(position_within_two)=[];
   Unique_interaction_counter=Unique_interaction_counter+1;
   
   %If center not within two fractions
   elseif nnz(position_within_two)==0
   interaction_final.unique_interactions(Unique_interaction_counter,1)= unique_interaction(interaction_counter2A);
   %check size of replicate_numbers
   interaction_final.replicate_numbers{Unique_interaction_counter,row_counter}=...
   cellfun(@str2num,binary_interaction_list(location_interaction_pairs(counter)+1,1));

   %copy the values for the first unique interaction to the list
   interaction_final.proteinA(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(1)+1,2);
   interaction_final.proteinB(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(1)+1,3);
   interaction_final.CenterA(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(1)+1,4);
   interaction_final.CenterB(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(1)+1,5);
   interaction_final.DeltaHeight(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(1)+1,6);
   interaction_final.DeltaCenter(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(1)+1,7);
   interaction_final.DeltaWidth(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(1)+1,8);
   interaction_final.DeltaEucDist(Unique_interaction_counter,row_counter)= binary_interaction_list(location_interaction_pairs(1)+1,9);   
   interaction_final.proteinInCorum(Unique_interaction_counter,1)= binary_interaction_list(location_interaction_pairs(1)+1,10);
   interaction_final.interactionInCorum(Unique_interaction_counter,1)= binary_interaction_list(location_interaction_pairs(1)+1,11);
   
   %Add one to counters and remove location already checked
   row_counter=row_counter+1; 
   location_interaction_pairs(1)=[];
   Unique_interaction_counter=Unique_interaction_counter+1;
   end
   end
end


%% True interactions Observation across replicate
%Count how many times a true positive interactions was observed
Total_unique_interactions=length(interaction_final.unique_interactions);
number_observation_pre_interaction=zeros(Total_unique_interactions,1);
number_unique_interaction=zeros(Total_unique_interactions,1);
%Record the number of times a interaction was detected in irrespective of replicate number
interactions_prep_replicate=zeros(Total_unique_interactions,replicate_counter);
for counter=1:Total_unique_interactions

%convert cell array to double
testA=[interaction_final.replicate_numbers{counter,:}];
%Save values
number_observation_pre_interaction(counter,1)=nnz(testA);

%Determine if multiple gaussians within the same replicate were indentified
unique_testA=unique(testA);
number_of_unique=length(unique_testA);
number_unique_interaction(counter,1)=number_of_unique;

%Which replicate is the interaction in
for replicate_counter=1:(number_of_replicates*number_of_channels)
interactions_prep_replicate(counter,replicate_counter)=nnz(find(testA(:)==replicate_counter));
end

end

%Copy number of unique_interaction across replicate to interaction array
interaction_final.times_observed_across_reps=number_unique_interaction;


%Determine how many times an interaction was observed
for replicate_counter=1:max(number_observation_pre_interaction)
number_observation_pre_replicate(replicate_counter)=nnz(find(number_observation_pre_interaction(:)==replicate_counter));
end

%Determine how many times an interaction was observed
for replicate_counter=1:max(number_unique_interaction)
number_unique_observation_pre_replicate(replicate_counter)=nnz(find(number_unique_interaction(:)==replicate_counter));
end

%total number of interactions of unique replicates
total_observation=sum(number_observation_pre_interaction);
%total number of interactions of unique replicates
total_unique_observation=sum(number_unique_interaction);


%% create unique identifier for false interactions
interaction_FN_pairs=strcat(FN_binary_interaction_list(:,2),'_',FN_binary_interaction_list(:,3));
unique_FN_interaction = unique(interaction_FN_pairs(2:end));
length_unique_FN_inter=length(unique_FN_interaction);

%Determine the number of unique false negatives
interactions_in_corum_FN=FN_binary_interaction_list(2:end,11);
interaction_as_number=cellfun(@str2double,interactions_in_corum_FN); 
%interaction_as_number=str2double(interactions_in_corum_FN);

%Create array to place values in
interaction_FN_incorum_pairs=cell(length(interactions_in_corum_FN),1);

parfor interaction_counter=1:length(interactions_in_corum_FN)
    if interaction_as_number(interaction_counter)==1
    interaction_FN_incorum_pairs(interaction_counter)=interaction_FN_pairs(interaction_counter+1);
    end
end

%remove empty cells
interaction_FN_incorum_pairs(cellfun('isempty',interaction_FN_incorum_pairs)) = [];

%Determine the unique interaction and count them
unique_FN_interaction_in_corum = unique(interaction_FN_incorum_pairs);
length_unique_FN_incorum=length(unique_FN_interaction_in_corum);

%clear varibles to speed up processing
clearvars interaction_FN_incorum_pairs;
clearvars interaction_FN_pairs;
clearvars FN_binary_interaction_list;

%% Caliculate TP, FP, FN and TN
%Caliculate values from true interaction data set 
True_proteinInCorum=sum([str2double(interaction_final.proteinInCorum)]);
True_interactionInCorum=sum([str2double(interaction_final.interactionInCorum)]);

%Caliculate values from true interaction data set 
False_proteinInCorum=length_unique_FN_inter;
False_interactionInCorum=length_unique_FN_incorum;

%caliculate TP, FP, FN and TN
final_TN=False_proteinInCorum-False_interactionInCorum;
final_FP=True_proteinInCorum-True_interactionInCorum;
final_FN=False_interactionInCorum;
final_TP=True_interactionInCorum;

% final results
final_Recall = final_TP/(final_TP+final_TN)
final_Precision = final_TP/(final_TP+final_FP)
final_TPR = final_TP/(final_TP+final_FN)
final_FPR =final_FP/(final_FP+final_TN)


%Format data to be written out
%Crete list of centerA, note ensure strjoin function is avalible 
interaction_final.centerA_formated=cell(Total_unique_interactions,1);

for format_loop=1:Total_unique_interactions
    centerA_2bformated1=interaction_final.CenterA(format_loop,:);
    centerA_2bformated1(cellfun('isempty',centerA_2bformated1)) = [];
    length_centerA=length(centerA_2bformated1);
    
    %Test if the array is longer then one entry
    if length_centerA<2
    interaction_final.centerA_formated(format_loop)=centerA_2bformated1(1);
    elseif length_centerA>=2
    interaction_final.centerA_formated(format_loop)= {strjoin(centerA_2bformated1,' ; ')}; 
    end
end

%Crete list of centerB, note ensure strjoin function is avalible 
interaction_final.centerB_formated=cell(Total_unique_interactions,1);

for format_loop=1:Total_unique_interactions
    centerB_2bformated1=interaction_final.CenterB(format_loop,:);
    centerB_2bformated1(cellfun('isempty',centerB_2bformated1)) = [];
    length_centerB=length(centerB_2bformated1);
    
    %Test if the array is longer then one entry
    if length_centerB<2
    interaction_final.centerB_formated(format_loop)=centerB_2bformated1(1);
    elseif length_centerB>=2
    interaction_final.centerB_formated(format_loop)= {strjoin(centerB_2bformated1,' ; ')}; 
    end
end

%Crete list of DeltaEuc, note ensure strjoin function is avalible 
interaction_final.DeltaEuc_formated=cell(Total_unique_interactions,1);

for format_loop=1:Total_unique_interactions
    DeltaEuc_2bformated1=interaction_final.DeltaEucDist(format_loop,:);
    DeltaEuc_2bformated1(cellfun('isempty',DeltaEuc_2bformated1)) = [];
    length_DeltaEuc=length(DeltaEuc_2bformated1);
    
    %Test if the array is longer then one entry
    if length_DeltaEuc<2
    interaction_final.DeltaEuc_formated(format_loop)=DeltaEuc_2bformated1(1);
    elseif length_DeltaEuc>=2
    interaction_final.DeltaEuc_formated(format_loop)= {strjoin(DeltaEuc_2bformated1,' ; ')}; 
    end
end

%Crete list of DeltaCenter, note ensure strjoin function is avalible 
interaction_final.DeltaCenter_formated=cell(Total_unique_interactions,1);

for format_loop=1:Total_unique_interactions
    DeltaCenter_2bformated1=interaction_final.DeltaCenter(format_loop,:);
    DeltaCenter_2bformated1(cellfun('isempty',DeltaCenter_2bformated1)) = [];
    length_DeltaCenter=length(DeltaCenter_2bformated1);
    
    %Test if the array is longer then one entry
    if length_DeltaCenter<2
    interaction_final.DeltaCenter_formated(format_loop)=DeltaCenter_2bformated1(1);
    elseif length_DeltaCenter>=2
    interaction_final.DeltaCenter_formated(format_loop)= {strjoin(DeltaCenter_2bformated1,' ; ')}; 
    end
end

%Crete list of DeltaHeight, note ensure strjoin function is avalible 
interaction_final.DeltaHeight_formated=cell(Total_unique_interactions,1);

for format_loop=1:Total_unique_interactions
    DeltaHeight_2bformated1=interaction_final.DeltaHeight(format_loop,:);
    DeltaHeight_2bformated1(cellfun('isempty',DeltaHeight_2bformated1)) = [];
    length_DeltaHeight=length(DeltaHeight_2bformated1);
    
    %Test if the array is longer then one entry
    if length_DeltaHeight<2
    interaction_final.DeltaHeight_formated(format_loop)=DeltaHeight_2bformated1(1);
    elseif length_DeltaHeight>=2
    interaction_final.DeltaHeight_formated(format_loop)= {strjoin(DeltaHeight_2bformated1,' ; ')}; 
    end
end

%Crete list of Deltawidth, note ensure strjoin function is avalible 
interaction_final.Deltawidth_formated=cell(Total_unique_interactions,1);

for format_loop=1:Total_unique_interactions
    Deltawidth_2bformated1=interaction_final.DeltaWidth(format_loop,:);
    Deltawidth_2bformated1(cellfun('isempty',Deltawidth_2bformated1)) = [];
    length_Deltawidth=length(Deltawidth_2bformated1);
    
    %Test if the array is longer then one entry
    if length_Deltawidth<2
    interaction_final.Deltawidth_formated(format_loop)=Deltawidth_2bformated1(1);
    elseif length_Deltawidth>=2
    interaction_final.Deltawidth_formated(format_loop)= {strjoin(Deltawidth_2bformated1,' ; ')}; 
    end
end

%Crete list of replicates, note ensure strjoin function is avalible 
interaction_final.Replicates_formated=cell(Total_unique_interactions,1);

for format_loop=1:Total_unique_interactions
    Replicates_2bformated1=interaction_final.replicate_numbers(format_loop,:);
    Replicates_2bformated1(cellfun('isempty',Replicates_2bformated1)) = [];
    length_Replicates=length(Replicates_2bformated1);
    
    %Test if the array is longer then one entry
    if length_Replicates<2
    interaction_final.Replicates_formated{format_loop}=mat2str(Replicates_2bformated1{1});
    elseif length_Replicates>=2
    interaction_final.Replicates_formated{format_loop}= strjoin([Replicates_2bformated1{:}],' ; '); 
    end
end


%create figure of precision replicates over replicates
%create array to fill with number of interaction with both proteins in
%corum and the number of interactions for these proteins in corum
Precision_array=zeros((number_of_replicates*number_of_channels),3);
Precision_redundant=zeros((number_of_replicates*number_of_channels),3);

%Determine the precision for interactions observed across replicates
for precision_counter=1:(number_of_replicates*number_of_channels)
   %Find position of interaction observed within the required number of
   %replicates
   position_to_sum=find(interaction_final.times_observed_across_reps(:)==precision_counter);

   %Count the number values for non redundant list
   Precision_array(precision_counter,1)=sum(cellfun(@str2double,interaction_final.proteinInCorum(position_to_sum)));
   Precision_array(precision_counter,2)=sum(cellfun(@str2double,interaction_final.interactionInCorum(position_to_sum)));
   Precision_array(precision_counter,3)=length(position_to_sum)
   %Count the number values for redundant list
   Precision_redundant(precision_counter,1)=(sum(cellfun(@str2double,interaction_final.proteinInCorum(position_to_sum)))*precision_counter);
   Precision_redundant(precision_counter,2)=(sum(cellfun(@str2double,interaction_final.interactionInCorum(position_to_sum)))*precision_counter);
   Precision_redundant(precision_counter,3)=(length(position_to_sum)*precision_counter);
end

%Generation precision as a precentage for non redundant interactions
Precent_precision_R=zeros((number_of_replicates*number_of_channels),1);
for precision_counter=1:(number_of_replicates*number_of_channels)
Precent_precision_R(precision_counter)=Precision_array(precision_counter,2)/(Precision_array(precision_counter,1));
end

%Generate precision as a precentage for redundant interactions
Precent_precision_NR=sum(Precision_redundant(:,2))/sum(Precision_redundant(:,1));

%% Write out True interaction and summary lists for precision.
%Generate headers to write out
size_interaction_replicates=size(interactions_prep_replicate);

%Create header for the list of replicate and interactions observed
for column_counter=1:size_interaction_replicates(2)
    Column_name1{column_counter}=strcat('Interaction Replicate_',mat2str(column_counter));
end

%Create header
Column_header1=repmat('%s,', 1, size_interaction_replicates(2)+7);

cd('Combined results');
Interaction_pre_rep_name=strcat('Interactions_across_replicate_',mat2str(Precision_values(precision_write_out_counter)),'pc.csv');
Inter_pre_rep = fopen(Interaction_pre_rep_name,'w');
fprintf (Inter_pre_rep,[Column_header1, '\n'],... %header for output
'Protein interaction','Center A','Center B',...
'Total number of observation of interaction',...
'Total number of replicates interaction observed in',Column_name1{:},...
'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for write_out_loop2=1:Total_unique_interactions   
fprintf(Inter_pre_rep,'%s,',interaction_final.unique_interactions{write_out_loop2,1});
fprintf(Inter_pre_rep,'%6.3f,',str2double(interaction_final.CenterA{write_out_loop2,1}));
fprintf(Inter_pre_rep,'%6.3f,',str2double(interaction_final.CenterB{write_out_loop2,1}));
fprintf(Inter_pre_rep,'%6f,',number_observation_pre_interaction(write_out_loop2,1));
fprintf(Inter_pre_rep,'%6f,', number_unique_interaction(write_out_loop2,1));
fprintf(Inter_pre_rep,'%6f,', interactions_prep_replicate(write_out_loop2,1:end));
fprintf(Inter_pre_rep,'%6f,',str2double(interaction_final.proteinInCorum{write_out_loop2,1}));
fprintf(Inter_pre_rep,'%6f,',str2double(interaction_final.interactionInCorum{write_out_loop2,1}));
fprintf(Inter_pre_rep,',\n');
end       
fclose(Inter_pre_rep);

Final_Results_name=strcat('Summary_Results_',mat2str(Precision_values(precision_write_out_counter)),'pc_replicate.csv');
fid_final= fopen(Final_Results_name,'wt'); % create the summary file of the interaction output 
fprintf (fid_final,'%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',... %header for OutputGaus output
                'Recall (non redundant)', 'Precision (non redundant)', 'TPR (non redundant)', 'FPR (non redundant)',...
                'Precision (redundant)','Total number of interactions (redundant across replicate)',...
                'Unique number of interactions (Nonredundant within replicate)',...
                'Total Unique interactions (Nonredundant across replicate)',...
                'Number non-redundant interactions' ); %Write Header
fprintf (fid_final,'%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,\n',...
                    final_Recall, final_Precision, final_TPR,...
                    final_FPR, Precent_precision_NR, total_observation,...
                    total_unique_observation, Total_unique_interactions, length_unique_inter);
fprintf (fid_final,',\n');             
fclose(fid_final);

Final_list_Interactionsname=strcat('Final_Interactions_list_',mat2str(Precision_values(precision_write_out_counter)),'_precision.csv');
fid_final_1 = fopen(Final_list_Interactionsname,'w');
fprintf (fid_final_1,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
        'Unique interactions','Protein A','Protein B','Center A',...
        'Center B','Replicate','Delta Height', 'Delta Center','Delta Width',...
        'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for ix=1:Total_unique_interactions    
fprintf(fid_final_1,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
interaction_final.unique_interactions{ix,1},...
interaction_final.proteinA {ix,1},...
interaction_final.proteinB {ix,1},...
interaction_final.centerA_formated{ix,1},...
interaction_final.centerB_formated{ix,1},...
interaction_final.Replicates_formated{ix,1},...
interaction_final.DeltaHeight_formated{ix,1},...
interaction_final.DeltaCenter_formated {ix,1},...
interaction_final.Deltawidth_formated{ix,1},...
interaction_final.DeltaEuc_formated{ix,1},...
interaction_final.proteinInCorum{ix,1},...
interaction_final.interactionInCorum{ix,1});
end
fclose(fid_final_1);

Final_Results_name=strcat('Global_Precision_across_replicates_',mat2str(Precision_values(precision_write_out_counter)),'pc.csv');
fid_final2= fopen(Final_Results_name,'wt'); % create the summary file of the interaction output 
fprintf (fid_final2,'%s,%s,%s,%s,%s,\n',... %header for OutputGaus output
                'Number of replicates interaction observed in',...
                'Total interactions number of interactions (non redundent interaction)',...
                'Interaction in corum (TP non redundent interaction)',...
                'Proteins in corum (TP + FP non redundent interaction)',...
                'Precision % (non redundent interaction)'); %Write Header
for precision_counter=1:(number_of_replicates*number_of_channels)
fprintf (fid_final2,'%6f,%6f,%6f,%6f,%6f,\n',...
                    precision_counter, Precision_array(precision_counter,3),...
                    Precision_array(precision_counter,2),...
                    Precision_array(precision_counter,1),...
                    Precent_precision_R(precision_counter));
end            
fclose(fid_final2);


%Plot Precision calculations
Interaction_not_in_corum=(Precision_array(:,3)-Precision_array(:,1));
Interaction_in_corum_not_detected=(Precision_array(:,1)-Precision_array(:,2));

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];

figure
f1=subplot(2,1,1)
f1_figure=bar(1:(number_of_replicates*number_of_channels), [Interaction_in_corum_not_detected(:,1) (Precision_array(:,2)) Interaction_not_in_corum(:,1)], 0.6, 'stack');

for k=1:3
set(f1_figure(k),'facecolor', myC(k,:), 'EdgeColor', 'k' );
end
legend('Proteins in corum (FP)', 'Interaction in corum (TP)',  'Proteins/Interactions not in corum');
ylim([0,(Interaction_in_corum_not_detected(1,1)+Precision_array(1,2)+Interaction_not_in_corum(1,1))*1.1]);
title('Interactions observed across isotoplogue channels','FontSize', 12);
ylabel('Number of interactions','FontSize', 8);
xlabel('isotoplogue channels','FontSize', 8);

f1=subplot(2,1,2)
f2_figure=bar(1:(number_of_replicates*number_of_channels), [Interaction_in_corum_not_detected(:,1) (Precision_array(:,2)) Interaction_not_in_corum(:,1)], 0.6, 'stack');

for k=1:3
set(f2_figure(k),'facecolor', myC(k,:), 'EdgeColor', 'k' );
end
legend('Proteins in corum (FP)', 'Interaction in corum (TP)',  'Proteins/Interactions not in corum');
ylim([0,(Interaction_in_corum_not_detected(3,1)+Precision_array(3,2)+Interaction_not_in_corum(3,1))*1.2]);
title('Interactions observed across isotoplogue channels (Zoom)','FontSize', 12);
ylabel('Number of interactions','FontSize', 8);
xlabel('isotoplogue channels','FontSize', 8);

% Save figure
Final_Interaction_figure=strcat('Observed_interactions_across_replicates_',mat2str(Precision_values(precision_write_out_counter)),'_precision.png');
saveas(f1, Final_Interaction_figure);


%Interaction observed only under treatment
treatment_replicates=[4,5,6];
untreatment_replicates=[1,2,3];

Observed_treatment=zeros(Total_unique_interactions,1);
Observed_untreatment=zeros(Total_unique_interactions,1);
treatment_specific_interaction=zeros(Total_unique_interactions,1);
untreatment_specific_interaction=zeros(Total_unique_interactions,1);

%Determine If there are condition specific interactions detected
for treatment_counter=1:Total_unique_interactions
    
   Observed_treatment(treatment_counter,1)=nnz(interactions_prep_replicate(treatment_counter,treatment_replicates));
   Observed_untreatment(treatment_counter,1)=nnz(interactions_prep_replicate(treatment_counter,untreatment_replicates));
   
   %Determine if there are treatment specific interactions detected
   if Observed_treatment(treatment_counter,1)>=2 & Observed_untreatment(treatment_counter,1)==0
      treatment_specific_interaction(treatment_counter,1)=1;
   end
   
   %Determine if there are untreatment specific interactions detected
   if Observed_untreatment(treatment_counter,1)>=2 & Observed_treatment(treatment_counter,1)==0
      untreatment_specific_interaction(treatment_counter,1)=1;
   end   
   
end

%Copy treatment specific interaction to list to write out
treatment_specific.unique_interactions=cell(0,0);
treatment_specific.proteinA=cell(0,0);
treatment_specific.proteinB=cell(0,0);
treatment_specific.centerA_formated=cell(0,0);
treatment_specific.centerB_formated=cell(0,0);
treatment_specific.Replicates_formated=cell(0,0);
treatment_specific.DeltaHeight_formated=cell(0,0);
treatment_specific.DeltaCenter_formated=cell(0,0);
treatment_specific.Deltawidth_formated=cell(0,0);
treatment_specific.DeltaEuc_formated=cell(0,0);
treatment_specific.proteinInCorum=cell(0,0);
treatment_specific.interactionInCorum=cell(0,0);

%Create counter
treatment_row_counter=1;
for treatment_counter=1:Total_unique_interactions
    if  treatment_specific_interaction(treatment_counter,1)==1;
    treatment_specific.unique_interactions(treatment_row_counter,1)=interaction_final.unique_interactions(treatment_counter,1);
    treatment_specific.proteinA(treatment_row_counter,1)=interaction_final.proteinA(treatment_counter,1);
    treatment_specific.proteinB(treatment_row_counter,1)=interaction_final.proteinB(treatment_counter,1);
    treatment_specific.centerA_formated(treatment_row_counter,1)=interaction_final.centerA_formated(treatment_counter,1);
    treatment_specific.centerB_formated(treatment_row_counter,1)=interaction_final.centerB_formated(treatment_counter,1);
    treatment_specific.Replicates_formated(treatment_row_counter,1)=interaction_final.Replicates_formated(treatment_counter,1);
    treatment_specific.DeltaHeight_formated(treatment_row_counter,1)=interaction_final.DeltaHeight_formated(treatment_counter,1);
    treatment_specific.DeltaCenter_formated(treatment_row_counter,1)=interaction_final.DeltaCenter_formated(treatment_counter,1);
    treatment_specific.Deltawidth_formated(treatment_row_counter,1)=interaction_final.Deltawidth_formated(treatment_counter,1);
    treatment_specific.DeltaEuc_formated(treatment_row_counter,1)=interaction_final.DeltaEuc_formated(treatment_counter,1);
    treatment_specific.proteinInCorum(treatment_row_counter,1)=interaction_final.proteinInCorum(treatment_counter,1);
    treatment_specific.interactionInCorum(treatment_row_counter,1)=interaction_final.interactionInCorum(treatment_counter,1);
    treatment_row_counter=1+treatment_row_counter;
    end
end

%Determine the number of interactions
treament_specific_interactionsnum=length(treatment_specific.unique_interactions);

%Write out treatment specific interactions
Final_list_Interactionsname=strcat('Final_Treatment_specific_interactions_list_',mat2str(Precision_values(precision_write_out_counter)),'_precision.csv');
fid_final_3 = fopen(Final_list_Interactionsname,'w');
fprintf (fid_final_3,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
        'Unique interactions','Protein A','Protein B','Center A',...
        'Center B','Replicate','Delta Height', 'Delta Center','Delta Width',...
        'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for xii=1:treament_specific_interactionsnum    
fprintf(fid_final_3,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
treatment_specific.unique_interactions{xii,1},...
treatment_specific.proteinA {xii,1},...
treatment_specific.proteinB {xii,1},...
treatment_specific.centerA_formated{xii,1},...
treatment_specific.centerB_formated{xii,1},...
treatment_specific.Replicates_formated{xii,1},...
treatment_specific.DeltaHeight_formated{xii,1},...
treatment_specific.DeltaCenter_formated {xii,1},...
treatment_specific.Deltawidth_formated{xii,1},...
treatment_specific.DeltaEuc_formated{xii,1},...
treatment_specific.proteinInCorum{xii,1},...
treatment_specific.interactionInCorum{xii,1});
end
fclose(fid_final_3);    

%Copy treatment specific interaction to list to write out
untreatment_specific.unique_interactions=cell(0,0);
untreatment_specific.proteinA=cell(0,0);
untreatment_specific.proteinB=cell(0,0);
untreatment_specific.centerA_formated=cell(0,0);
untreatment_specific.centerB_formated=cell(0,0);
untreatment_specific.Replicates_formated=cell(0,0);
untreatment_specific.DeltaHeight_formated=cell(0,0);
untreatment_specific.DeltaCenter_formated=cell(0,0);
untreatment_specific.Deltawidth_formated=cell(0,0);
untreatment_specific.DeltaEuc_formated=cell(0,0);
untreatment_specific.proteinInCorum=cell(0,0);
untreatment_specific.interactionInCorum=cell(0,0);

%Create counter
untreatment_row_counter=1;
for treatment_counter=1:Total_unique_interactions
    if  untreatment_specific_interaction(treatment_counter,1)==1;
    untreatment_specific.unique_interactions(untreatment_row_counter,1)=interaction_final.unique_interactions(treatment_counter,1);
    untreatment_specific.proteinA(untreatment_row_counter,1)=interaction_final.proteinA(treatment_counter,1);
    untreatment_specific.proteinB(untreatment_row_counter,1)=interaction_final.proteinB(treatment_counter,1);
    untreatment_specific.centerA_formated(untreatment_row_counter,1)=interaction_final.centerA_formated(treatment_counter,1);
    untreatment_specific.centerB_formated(untreatment_row_counter,1)=interaction_final.centerB_formated(treatment_counter,1);
    untreatment_specific.Replicates_formated(untreatment_row_counter,1)=interaction_final.Replicates_formated(treatment_counter,1);
    untreatment_specific.DeltaHeight_formated(untreatment_row_counter,1)=interaction_final.DeltaHeight_formated(treatment_counter,1);
    untreatment_specific.DeltaCenter_formated(untreatment_row_counter,1)=interaction_final.DeltaCenter_formated(treatment_counter,1);
    untreatment_specific.Deltawidth_formated(untreatment_row_counter,1)=interaction_final.Deltawidth_formated(treatment_counter,1);
    untreatment_specific.DeltaEuc_formated(untreatment_row_counter,1)=interaction_final.DeltaEuc_formated(treatment_counter,1);
    untreatment_specific.proteinInCorum(untreatment_row_counter,1)=interaction_final.proteinInCorum(treatment_counter,1);
    untreatment_specific.interactionInCorum(untreatment_row_counter,1)=interaction_final.interactionInCorum(treatment_counter,1);
    untreatment_row_counter=1+untreatment_row_counter;
    end
end

%Determine the number of interactions
untreament_specific_interactionsnum=length(untreatment_specific.unique_interactions);

%Write out treatment specific interactions
Final_list_Interactionsname=strcat('Final_Untreatment_specific_interactions_list_',mat2str(Precision_values(precision_write_out_counter)),'_precision.csv');
fid_final_4 = fopen(Final_list_Interactionsname,'w');
fprintf (fid_final_4,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
        'Unique interactions','Protein A','Protein B','Center A',...
        'Center B','Replicate','Delta Height', 'Delta Center','Delta Width',...
        'Delta EucDist', 'Both proteins in Corum', 'Interaction in Corum'); %Write Header
for xiii=1:untreament_specific_interactionsnum    
fprintf(fid_final_4,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,\n',...
untreatment_specific.unique_interactions{xiii,1},...
untreatment_specific.proteinA {xiii,1},...
untreatment_specific.proteinB {xiii,1},...
untreatment_specific.centerA_formated{xiii,1},...
untreatment_specific.centerB_formated{xiii,1},...
untreatment_specific.Replicates_formated{xiii,1},...
untreatment_specific.DeltaHeight_formated{xiii,1},...
untreatment_specific.DeltaCenter_formated {xiii,1},...
untreatment_specific.Deltawidth_formated{xiii,1},...
untreatment_specific.DeltaEuc_formated{xiii,1},...
untreatment_specific.proteinInCorum{xiii,1},...
untreatment_specific.interactionInCorum{xiii,1});
end
fclose(fid_final_4);  
cd('../');
end

fclose('all');
