clear all

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

%% Initialize variables.
filename = 'Single_peptide_unique_Scaffold_peptide_list_reverse_contaminates_removed_With_Quant.txt';

%% Open the text file.
fileID = fopen(filename,'r');

%% Read columns of data according to format string.
dataArray = textscan(fileID, '%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%[^\n\r]',...
            'Delimiter', '\t', 'HeaderLines' ,1, 'ReturnOnError', false);

%% Close the text file.
fclose(fileID);


%% Allocate imported array to column variable names
NTAILS.Experimentname = dataArray{:, 1};
NTAILS.Biologicalsamplecategory = dataArray{:, 2};
NTAILS.Biologicalsamplename = dataArray{:, 3};
NTAILS.MSMSsamplename = dataArray{:, 4};
NTAILS.Proteingroup = dataArray{:, 5};
NTAILS.Proteinaccessionnumber = dataArray{:, 6};
NTAILS.UniprotTerms = dataArray{:, 7};
NTAILS.Proteinname = dataArray{:, 8};
NTAILS.Proteinidentificationprobability = dataArray{:, 9};
NTAILS.Proteinpercentageoftotalspectra = dataArray{:, 10};
NTAILS.ProteinmolecularweightAMU = dataArray{:, 11};
NTAILS.Numberofuniquepeptides = dataArray{:, 12};
NTAILS.Numberofuniquespectra = dataArray{:, 13};
NTAILS.Numberoftotalspectra = dataArray{:, 14};
NTAILS.Spectrumname = dataArray{:, 15};
NTAILS.Peptidesequence = dataArray{:, 16};
NTAILS.PeptideSequencewithmodifications = dataArray{:, 17};
NTAILS.Previousaminoacid = dataArray{:, 18};
NTAILS.Nextaminoacid = dataArray{:, 19};
NTAILS.Peptideidentificationprobability = dataArray{:, 20};
NTAILS.MascotIonscore = dataArray{:, 21};
NTAILS.MascotIdentityscore = dataArray{:, 22};
NTAILS.MascotDeltaIonscore = dataArray{:, 23};
NTAILS.Modificationsidentifiedbyspectrum = dataArray{:, 24};
NTAILS.ActualpeptidemassAMU = dataArray{:, 25};
NTAILS.Spectrumcharge = dataArray{:, 26};
NTAILS.ActualminuscalculatedpeptidemassAMU = dataArray{:, 27};
NTAILS.ActualminuscalculatedpeptidemassPPM = dataArray{:, 28};
NTAILS.Peptidestartindex = dataArray{:, 29};
NTAILS.Peptidestopindex = dataArray{:, 30};
NTAILS.Replicate1HvsL = dataArray{:, 31};
NTAILS.Replicate1MvsL = dataArray{:, 32};
NTAILS.Replicate1MvsH = dataArray{:, 33};
NTAILS.Replicate2HvsL = dataArray{:, 34};
NTAILS.Replicate2MvsL = dataArray{:, 35};
NTAILS.Replicate2MvsH = dataArray{:, 36};
NTAILS.Replicate3HvsL = dataArray{:, 37};
NTAILS.Replicate3MvsL = dataArray{:, 38};
NTAILS.Replicate3MvsH = dataArray{:, 39};
NTAILS.Replicate1NormalizedHvsL = dataArray{:, 40};
NTAILS.Replicate1NormalizedMvsL = dataArray{:, 41};
NTAILS.Replicate1NormalizedMvsH = dataArray{:, 42};
NTAILS.Replicate2NormalizedHvsL = dataArray{:, 43};
NTAILS.Replicate2NormalizedMvsL = dataArray{:, 44};
NTAILS.Replicate2NormalizedMvsH = dataArray{:, 45};
NTAILS.Replicate3NormalizedHvsL = dataArray{:, 46};
NTAILS.Replicate3NormalizedMvsL = dataArray{:, 47};
NTAILS.Replicate3NormalizedMvsH = dataArray{:, 48};
NTAILS.AverageNormalizedHvsL = dataArray{:, 49};
NTAILS.AverageNormalizedMvsL = dataArray{:, 50};
NTAILS.ttestpvalueHvsL = dataArray{:, 51};
NTAILS.ttestpvalueMvsL = dataArray{:, 52};
NTAILS.BenjaminiHochbergcorrectedpvalueHvsLttest = dataArray{:, 53};
NTAILS.BenjaminiHochbergcorrectedpvalueMvsLttest = dataArray{:, 54};
NTAILS.BenjaminiHochbergcorrectedpvalueHvsLMWW = dataArray{:, 55};
NTAILS.BenjaminiHochbergcorrectedpvalueMvsLMWW = dataArray{:, 56};


%% Initialize variables.
filename = 'Mapped_terms_Ntails.txt';
%% Open the text file.
fileID = fopen(filename,'r');

%% Read columns of data according to format string.
dataArray = textscan(fileID, '%s%s%s%s%s%s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%[^\n\r]', 'Delimiter', '\t',  'ReturnOnError', false);

%% Close the text file.
fclose(fileID);

%% Allocate imported array to column variable names
Mapped_terms.yourlistM2015071097E05OJGPJ = dataArray{:, 1};
Mapped_terms.Entry = dataArray{:, 2};
Mapped_terms.Entryname = dataArray{:, 3};
Mapped_terms.Status = dataArray{:, 4};
Mapped_terms.Proteinnames = dataArray{:, 5};
Mapped_terms.Genenames = dataArray{:, 6};

%Map unitprot accession numbers to proteins
NTAILS.uniprot_accession_num=cell(length(NTAILS.Experimentname),1);

for entry_counter=1:length(NTAILS.Experimentname)
    %check if any of the protein entries matahc the tails data
    [Position_of_uniprot_accession_num]=ind2sub(size(Mapped_terms.Entryname), strmatch(NTAILS.UniprotTerms{entry_counter}, Mapped_terms.Entryname, 'exact'));
    
    %record results
    if ~isempty(Position_of_uniprot_accession_num)
    NTAILS.uniprot_accession_num(entry_counter)=Mapped_terms.Entry(Position_of_uniprot_accession_num);
    else
    NTAILS.uniprot_accession_num(entry_counter)={''};    
    end
end

%Determine number of unique protein in tails data
Unique_N_tails_uniprot=unique(NTAILS.uniprot_accession_num);

%Remove entries without match to uniprot
Unique_N_tails_uniprot(cellfun(@isempty, Unique_N_tails_uniprot))=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Import master Gaussian list
filename = 'Master_guassian_list_cyto.csv';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, '%s%s%f%s%f%f%f%f%f%f%f%f%[^\n\r]', 'Delimiter', ',', 'EmptyValue' ,NaN,'HeaderLines' ,1, 'ReturnOnError', false);

%% Close the text file.
fclose(fileID);

Master_Gaussian_cyto.Proteinname = dataArray{:, 1};
Master_Gaussian_cyto.Unique_identifierofbestgaussian = dataArray{:, 2};
Master_Gaussian_cyto.Replicateofbestgaussian = dataArray{:, 3};
Master_Gaussian_cyto.Channelofbestgaussian = dataArray{:, 4};
Master_Gaussian_cyto.Guassian_index_numberofbestgaussian = dataArray{:, 5};
Master_Gaussian_cyto.Centerofbestgaussian = dataArray{:, 6};
Master_Gaussian_cyto.Heightofbestgaussian = dataArray{:, 7};
Master_Gaussian_cyto.Widthofbestgaussian = dataArray{:, 8};
Master_Gaussian_cyto.SSEofbestgaussian = dataArray{:, 9};
Master_Gaussian_cyto.adjrsquareofbestgaussian = dataArray{:, 10};
Master_Gaussian_cyto.GaussianAreaofbestgaussian = dataArray{:, 11};
Master_Gaussian_cyto.ComplexSizeofbestgaussian = dataArray{:, 12};

number_unique_guassian=length(Master_Gaussian_cyto.Proteinname);

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

fid_annotation_cyto=fopen('Master_guassian_list_fold_comparsion_cytoplasmic.csv', 'rt');                                                                  
Import_annotation_cyto= textscan(fid_annotation_cyto, '%s\b\t', 'Delimiter','\n', 'treatAsEmpty',{'N/A'});

%unnest import file
Import_annotation_cyto=Import_annotation_cyto{1,1};

%count number of cells
No_of_cells=length(Import_annotation_cyto);

%split file based on ,
for counter=1:No_of_cells
seperated_values=regexp(Import_annotation_cyto{counter,1},',','split');
for counter2=1:31
Imported_seperated_cyto(counter, counter2)=seperated_values(counter2);
end
end

%Convert values to usible form
Cyto_normalized_values=zeros(No_of_cells-1,1);
for counter=2:No_of_cells
    
Cyto_normalized_values(counter-1)=str2num(Imported_seperated_cyto{counter,18});
end

Proteins_with_Gaussian=unique(Imported_seperated_cyto(2:end,1));
Number_proteins_with_gaussian=length(Proteins_with_Gaussian);

%% Import master Gaussian list
filename = 'Master_guassian_list_mito.csv';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, '%s%s%f%s%f%f%f%f%f%f%f%f%[^\n\r]', 'Delimiter', ',', 'EmptyValue' ,NaN,'HeaderLines' ,1, 'ReturnOnError', false);

%% Close the text file.
fclose(fileID);

Master_Gaussian_mito.Proteinname = dataArray{:, 1};
Master_Gaussian_mito.Unique_identifierofbestgaussian = dataArray{:, 2};
Master_Gaussian_mito.Replicateofbestgaussian = dataArray{:, 3};
Master_Gaussian_mito.Channelofbestgaussian = dataArray{:, 4};
Master_Gaussian_mito.Guassian_index_numberofbestgaussian = dataArray{:, 5};
Master_Gaussian_mito.Centerofbestgaussian = dataArray{:, 6};
Master_Gaussian_mito.Heightofbestgaussian = dataArray{:, 7};
Master_Gaussian_mito.Widthofbestgaussian = dataArray{:, 8};
Master_Gaussian_mito.SSEofbestgaussian = dataArray{:, 9};
Master_Gaussian_mito.adjrsquareofbestgaussian = dataArray{:, 10};
Master_Gaussian_mito.GaussianAreaofbestgaussian = dataArray{:, 11};
Master_Gaussian_mito.ComplexSizeofbestgaussian = dataArray{:, 12};

number_unique_guassian_mito=length(Master_Gaussian_mito.Proteinname);

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

%Import majority protein IDs from membrane search
[~, Protein_IDs_mito] = xlsread('MajorProteinIDs_mito.xlsx');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Open the text file.
fileID = fopen('Peptide_Summary_HvsM_table.csv','r');

%% Read columns of data according to format string.
% This call is based on the structure of the file used to generate this
% code. If an error occurs for a different file, try regenerating the code
% from the Import Tool.
dataArray = textscan(fileID, '%s%s%f%s%s%s%s%f%f%[^\n\r]', 'Delimiter', ',', 'EmptyValue' ,NaN,'HeaderLines' ,1, 'ReturnOnError', false);

%% Allocate imported array to column variable names
diversity_in_peptide.Protein = dataArray{:, 1};
diversity_in_peptide.Isotopicchannel = dataArray{:, 2};
diversity_in_peptide.Replicate = dataArray{:, 3};
diversity_in_peptide.Containspeptidewhichdifferbasedonstdevrawvaluesfrommaxquant = dataArray{:, 6};
diversity_in_peptide.Containspeptidewhichdifferbasedonstdevfixedvalues = dataArray{:, 7};
diversity_in_peptide.ShannondiversityindexH = dataArray{:, 8};
diversity_in_peptide.Evenness = dataArray{:, 9};

Unique_proteins_Diversity=unique(diversity_in_peptide.Protein);

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

%% Initialize variables.
filename = 'Perseus_enrichment_Protein_level_file.csv';

%% Open the text file.
fileID = fopen(filename,'r');
dataArray = textscan(fileID, '%q%q%q%q%q%q%[^\n\r]', 'Delimiter', ',', 'EmptyValue' ,NaN, 'ReturnOnError', false);

%% Close the text file.
fclose(fileID);

Perseus_enrichment_Protein_level.Proteinname = dataArray{:, 1};
Perseus_enrichment_Protein_level.IncreaseANDDecreasestworeplicates = dataArray{:, 2};
Perseus_enrichment_Protein_level.Increasetworeplicates = dataArray{:, 3};
Perseus_enrichment_Protein_level.Decreasetworeplicates = dataArray{:, 4};
Perseus_enrichment_Protein_level.ttest_pvaluebonferronicorrectio = dataArray{:, 5};
Perseus_enrichment_Protein_level.MWW_pvaluebonferronicorrecti1 = dataArray{:, 6};

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

%%%Import degrabase
[~, ~, degrabase] = xlsread('degrabase_Version_1_downloaded_30012014.xlsx','Sheet1');
degrabase(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),degrabase)) = {''};

%define dimensions of degrabase
degrabase_dimension=size(degrabase);

%create array to copy names to
degrbase_names=cell(degrabase_dimension(1)-1,1);

%Convert degrabase names to array
for counter1=2:degrabase_dimension(1)
name_to_copy=degrabase(counter1,1);
degrbase_names(counter1-1,1)=name_to_copy;
end

%create array to copy names to
degrbase_P1_position=cell(degrabase_dimension(1)-1,1);

%Convert degrabase names to array
for counter1=2:degrabase_dimension(1)
name_to_P1=degrabase(counter1,10);
degrbase_P1_position(counter1-1,1)=name_to_P1;
end

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


%Determine which proteins in the N-terminome show chnages
Signficant_termini_Persaus=cell(length(Unique_N_tails_uniprot),2);

for Nterminal=1:length(Unique_N_tails_uniprot)
 Signficant_terminiMWW=0;
 Signficant_termini_ttest=0;
 
 %Reset position value
 Position1=[];
 %Look up uniprot name
[Position1]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Unique_N_tails_uniprot{Nterminal}, NTAILS.uniprot_accession_num,'exact'));

if ~isempty(Position1)
    
    for unique_counter=1:length(Position1) 
    %check if any of the termini show signficant chnages
    if str2double(NTAILS.BenjaminiHochbergcorrectedpvalueMvsLMWW(Position1(unique_counter)))<0.05 |...
       str2double(NTAILS.BenjaminiHochbergcorrectedpvalueHvsLMWW(Position1(unique_counter)))<0.05 
       Signficant_terminiMWW=1;    
    end
    
    %check if any of the termini show signficant chnages
    if str2double(NTAILS.BenjaminiHochbergcorrectedpvalueMvsLttest(Position1(unique_counter)))<0.05 |...
       str2double(NTAILS.BenjaminiHochbergcorrectedpvalueHvsLttest(Position1(unique_counter)))<0.05 
       Signficant_termini_ttest=1;    
    end
    end
    
    if Signficant_terminiMWW==1
    Signficant_termini_Persaus(Nterminal,1)={'+'};
    else
    Signficant_termini_Persaus(Nterminal,1)={''};
    end
    
    if Signficant_termini_ttest==1
    Signficant_termini_Persaus(Nterminal,2)={'+'};
    else
    Signficant_termini_Persaus(Nterminal,2)={''};
    end    
    

end

end

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

%Check if in Degrabase at Protein level
IN_degrabase_Protein=cell(length(Unique_N_tails_uniprot),1);
Protein_counter=0;
for Nterminal=1:length(Unique_N_tails_uniprot)
    
 %Reset position value
 location_Protein_in_degabase=[];
 
    [location_Protein_in_degabase]=ind2sub(size(degrbase_names), strmatch(Unique_N_tails_uniprot{Nterminal},degrbase_names(:),'exact'));

    if ~isempty(location_Protein_in_degabase)
    IN_degrabase_Protein(Nterminal)={'+'};
    Protein_counter=Protein_counter+1;
    else
    IN_degrabase_Protein(Nterminal,2)={''};
    end
end

%Check if in Degrabase at cleavage level
IN_degrabase_site=cell(length(NTAILS.uniprot_accession_num),1);
Site_counter=0;

for Degrabase_cleavage_level=1:length(NTAILS.uniprot_accession_num)
%Reset observed value
Observed_value=0;

%Find position
[Position1]=ind2sub(size(degrabase(:,1)), strmatch(NTAILS.uniprot_accession_num{Degrabase_cleavage_level}, degrabase(:,1),'exact'));

for position_counter=1:length(Position1)
    if isequal(degrabase{Position1(position_counter),6},str2double(NTAILS.Peptidestartindex(Degrabase_cleavage_level))) & Observed_value==0
    IN_degrabase_site(Degrabase_cleavage_level)={'+'};
    Site_counter=Site_counter+1;
    Observed_value=1;
    end
end


end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Save Scaffold temp non-redunent list of peptides with Quant values
%Write out combined peptide list
fid = fopen('Non_redundent_unique_Scaffold_peptide_list_reverse_contaminates_removed_With_Quant_mapped.txt','w');
fprintf (fid,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n',... %header for output
'Experiment name','Biological sample category','Biological sample name','MS/MS sample name','Protein group',...
'Protein accession number','Uniprot Terms','Uniprot entry','Protein name','Protein identification probability','Protein percentage of total spectra',...
'Protein molecular weight (AMU)','Number of unique peptides','Number of unique spectra','Number of total spectra',...
'Spectrum name','Peptide sequence','Peptide Sequence(with modifications)','Previous amino acid','Next amino acid','Peptide identification probability',...
'Mascot Ion score','Mascot Identity score','Mascot Delta Ion score','Modifications identified by spectrum',...
'Actual peptide mass (AMU)','Spectrum charge','Actual minus calculated peptide mass (AMU)',...
'Actual minus calculated peptide mass (PPM)','Peptide start index','Peptide stop index', 'Replicate 1 HvsL',...
'Replicate 1 MvsL','Replicate 1 MvsH','Replicate 2 HvsL','Replicate 2 MvsL','Replicate 2 MvsH',...
'Replicate 3 HvsL','Replicate 3 MvsL','Replicate 3 MvsH','Replicate 1 Normalized HvsL',...
'Replicate 1 Normalized MvsL','Replicate 1 Normalized MvsH','Replicate 2 Normalized HvsL',...
'Replicate 2 Normalized MvsL','Replicate 2 Normalized MvsH',...
'Replicate 3 Normalized HvsL','Replicate 3 Normalized MvsL','Replicate 3 Normalized MvsH',...
'Average Normalized HvsL','Average Normalized MvsL','ttest p-value HvsL',....
'ttest p-value MvsL', 'Benjamini-Hochberg corrected p-value HvsL (ttest)',...
'Benjamini-Hochberg corrected p-value MvsL (ttest)','Benjamini-Hochberg corrected p-value HvsL (MWW)',...
'Benjamini-Hochberg corrected p-value MvsL (MWW)','Site known in Degrabase'); %Write Header
for write_out_loop2=1:length(NTAILS.Experimentname)
%Write out data          
fprintf(fid,'%s	',NTAILS.Experimentname{write_out_loop2});
fprintf(fid,'%s	',NTAILS.Biologicalsamplecategory{write_out_loop2});
fprintf(fid,'%s	',NTAILS.Biologicalsamplename{write_out_loop2});
fprintf(fid,'%s	',NTAILS.MSMSsamplename{write_out_loop2});
fprintf(fid,'%s	',NTAILS.Proteingroup{write_out_loop2});
fprintf(fid,'%s	',NTAILS.Proteinaccessionnumber{write_out_loop2});
fprintf(fid,'%s	',NTAILS.UniprotTerms{write_out_loop2});
fprintf(fid,'%s	',NTAILS.uniprot_accession_num{write_out_loop2});
fprintf(fid,'%s	',NTAILS.Proteinname{write_out_loop2});
fprintf(fid,'%f	',str2double(NTAILS.Proteinidentificationprobability{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Proteinpercentageoftotalspectra{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.ProteinmolecularweightAMU{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Numberofuniquepeptides{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Numberofuniquespectra{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Numberoftotalspectra{write_out_loop2}));
fprintf(fid,'%s	',NTAILS.Spectrumname{write_out_loop2});
fprintf(fid,'%s	',NTAILS.Peptidesequence{write_out_loop2});
fprintf(fid,'%s	',NTAILS.PeptideSequencewithmodifications{write_out_loop2});
fprintf(fid,'%s	',NTAILS.Previousaminoacid {write_out_loop2});
fprintf(fid,'%s	',NTAILS.Nextaminoacid{write_out_loop2});
fprintf(fid,'%f	',str2double(NTAILS.Peptideidentificationprobability{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.MascotIonscore{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.MascotIdentityscore{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.MascotDeltaIonscore{write_out_loop2}));
fprintf(fid,'%s	',NTAILS.Modificationsidentifiedbyspectrum{write_out_loop2});
fprintf(fid,'%f	',str2double(NTAILS.ActualpeptidemassAMU{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Spectrumcharge{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.ActualminuscalculatedpeptidemassAMU{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.ActualminuscalculatedpeptidemassPPM{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Peptidestartindex{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Peptidestopindex{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate1HvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate1MvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate1MvsH{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate2HvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate2MvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate2MvsH{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate3HvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate3MvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate3MvsH{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate1NormalizedHvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate1NormalizedMvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate1NormalizedMvsH{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate2NormalizedHvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate2NormalizedMvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate2NormalizedMvsH{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate3NormalizedHvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate3NormalizedMvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.Replicate3NormalizedMvsH{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.AverageNormalizedHvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.AverageNormalizedMvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.ttestpvalueHvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.ttestpvalueMvsL{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.BenjaminiHochbergcorrectedpvalueHvsLttest{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.BenjaminiHochbergcorrectedpvalueMvsLttest{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.BenjaminiHochbergcorrectedpvalueHvsLMWW{write_out_loop2}));
fprintf(fid,'%f	',str2double(NTAILS.BenjaminiHochbergcorrectedpvalueMvsLMWW{write_out_loop2}));
fprintf(fid,'%s	',IN_degrabase_site{write_out_loop2});
fprintf(fid,'\n');
end
fclose(fid);


%Write out table for enrichment analysis
Output_name=strcat('Persaus_Protein_level_Ntails.txt');

fid = fopen(Output_name,'w');
fprintf (fid,'%s\t%s\t%s\t%s\t\n',...
'Protein name',...
'N-tail change below adjusted_pvalue(Ben Hochberg) (ttest)',....
'N-tail change below adjusted_pvalue(Ben Hochberg) (MWW)',...
'Protein within Degrabase'); %Write Header
for counter_output=1:length(Unique_N_tails_uniprot)
fprintf(fid,'%s	',Unique_N_tails_uniprot{counter_output});
fprintf(fid,'%s	',Signficant_termini_Persaus{counter_output,2});
fprintf(fid,'%s	',Signficant_termini_Persaus{counter_output,1});
fprintf(fid,'%s	',IN_degrabase_Protein{counter_output,1});
fprintf(fid,'\n');
end
fclose(fid);

%Write out table for enrichment analysis
Output_name=strcat('Coverage_of_termini_and_protein_compared_to_degrabase.txt');

fid = fopen(Output_name,'w');
fprintf (fid,'%s\t%s\t%s\t%s\t\n',...
'Total number of observed Protein in tails dataset',...
'Total number of Termini in tails dataset',...
'Number of observed proteins in degrabase',....
'Nunmber of sites in degrabase'); %Write Header
fprintf(fid,'%f	',length(Unique_N_tails_uniprot));
fprintf(fid,'%f	',length(unique(NTAILS.Peptidesequence)));
fprintf(fid,'%f	',Protein_counter);
fprintf(fid,'%f	',Site_counter);
fprintf(fid,'\n');
fclose(fid);

%Determine if there is a correlation with proteins that change between the
%two type of MS data
%work out unique proteins in Diversity array
Unique_proteins_Diversity=unique(diversity_in_peptide.Protein);

Overlap_array_for_correlation=zeros(1,4);
Proteins_which_overlap=cell(1,2);
Number_of_terminin_per_proteins=zeros(1,1);
Signficant_termini=cell(1,2);
counter=1;


for Protein_counter1=1:length(Unique_proteins_Diversity)
    
    %Find protein in array 
    %check if any of the protein entries matahc the tails data
    [Position1]=ind2sub(size(diversity_in_peptide.Protein), strmatch(Unique_proteins_Diversity{Protein_counter1}, diversity_in_peptide.Protein,'exact'));
    Temp_value=nanmean(diversity_in_peptide.Evenness(Position1));
    Temp_value2=nanmean(diversity_in_peptide.ShannondiversityindexH(Position1));
    
    %Find protein in list of protein groups
    [PositionA, PositionB]=ind2sub(size(Protein_IDs), strmatch(Unique_proteins_Diversity{Protein_counter1}, Protein_IDs,'exact'));
    
    Checked_row_counter=0;
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs(PositionA,:)))
    
    %check if any of the protein entries matahc the tails data
    PositionC=[];
    [PositionC]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Protein_IDs{PositionA,terms_in_protein_group}, NTAILS.uniprot_accession_num));

    %test if
    if ~isempty(PositionC) & Checked_row_counter==0
    
    %determine the unique peptide
    A_1=unique(NTAILS.Peptidesequence(PositionC));
    B_1=NTAILS.Peptidesequence(PositionC);
    %Find position of unique terms
    [~,unique_temp]=ismember(A_1,B_1);
    
    Temp_value2_HvsL=NaN(1,1);
    Temp_value3_MvsL=NaN(1,1);
    Signficant_terminiMWW=NaN(1,1);
    Signficant_termini_ttest=NaN(1,1);
    
    for unique_counter=1:length(unique_temp)
    
    %Determine average peptide measurements
    Temp_value2_HvsL(unique_counter)=str2double(NTAILS.AverageNormalizedHvsL(PositionC(unique_temp(unique_counter))));
    Temp_value3_MvsL(unique_counter)=str2double(NTAILS.AverageNormalizedMvsL(PositionC(unique_temp(unique_counter))));
    
    %check if any of the termini show signficant chnages
    if str2double(NTAILS.BenjaminiHochbergcorrectedpvalueMvsLMWW(PositionC(unique_temp(unique_counter))))<0.05 |...
       str2double(NTAILS.BenjaminiHochbergcorrectedpvalueHvsLMWW(PositionC(unique_temp(unique_counter))))<0.05 
       Signficant_terminiMWW=1;    
    end
    
    %check if any of the termini show signficant chnages
    if str2double(NTAILS.BenjaminiHochbergcorrectedpvalueMvsLttest(PositionC(unique_temp(unique_counter))))<0.05 |...
       str2double(NTAILS.BenjaminiHochbergcorrectedpvalueHvsLttest(PositionC(unique_temp(unique_counter))))<0.05 
       Signficant_termini_ttest=1;    
    end
    
    end
    
    %Write out values
    Proteins_which_overlap(counter,1)=Unique_proteins_Diversity(Protein_counter1);
    Proteins_which_overlap(counter,2)=NTAILS.uniprot_accession_num(PositionC(1));
    Overlap_array_for_correlation(counter,1)=Temp_value;
    Overlap_array_for_correlation(counter,4)=Temp_value2;
    Overlap_array_for_correlation(counter,3)=nanmean(Temp_value2_HvsL);
    Overlap_array_for_correlation(counter,2)=nanmean(Temp_value3_MvsL);
    Number_of_terminin_per_proteins(counter)=length(unique_temp);
    
    if Signficant_terminiMWW==1
    Signficant_termini(counter,1)={'+'};
    else
    Signficant_termini(counter,1)={''};
    end
    
    if Signficant_termini_ttest==1
    Signficant_termini(counter,2)={'+'};
    else
    Signficant_termini(counter,2)={''};
    end    
    
    counter=counter+1;
    Checked_row_counter=1;
        
    end
    end
end
    

%determine the Number of unique uniprot terms which we observe gaussians for and look up in the N-tails
%dataset
Unique_proteins_cyto=unique(Master_Gaussian_cyto.Proteinname);
Overlap_array1=cell(1,2);
Overlap_counter1=1;
Not_overlap_counter=0;

for unique_counter=1:length(Unique_proteins_cyto)
    
    %look protein in major proein group file
    [row col]=ind2sub(size(Protein_IDs), strmatch(Unique_proteins_cyto{unique_counter}, Protein_IDs,'exact')); 
    
    Checked_row_counter=0;
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs(row,:)))
        
    %check if any of the protein entries matahc the tails data
    [Position]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Protein_IDs{row,terms_in_protein_group}, NTAILS.uniprot_accession_num, 'exact'));

    %test if
    if ~isempty(Position) & Checked_row_counter==0
        Overlap_array1(Overlap_counter1,1)=Unique_proteins_cyto(unique_counter);
        Overlap_array1(Overlap_counter1,2)=NTAILS.uniprot_accession_num(Position(1));
        Overlap_counter1=1+Overlap_counter1;  
        Checked_row_counter=1;
    else
        Not_overlap_counter=1+Not_overlap_counter;
    end
    end
    
end


%determine the Number of unique uniprot terms which we observe gaussians for and look up in the N-tails
%dataset
Overlap_array2=cell(1,2);
Overlap_counter2=1;

for unique_counter=2:length(Protein_IDs)
    Checked_row_counter=0;    
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs(unique_counter,:)))
    %check if any of the protein entries matahc the tails data
    [Position]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Protein_IDs{unique_counter,terms_in_protein_group}, NTAILS.uniprot_accession_num, 'exact'));

    %test if
    if ~isempty(Position) & Checked_row_counter==0
        Overlap_array2(Overlap_counter2,1)=Protein_IDs(unique_counter);
        Overlap_array2(Overlap_counter2,2)=NTAILS.uniprot_accession_num(Position(1));
        Overlap_counter2=1+Overlap_counter2;
        Checked_row_counter=1;
    end
     end
end

Overlap_array2_mito=cell(1,2);
Overlap_counter2_mito=1;

for unique_counter=2:length(Protein_IDs_mito)
    Checked_row_counter=0;    
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs_mito(unique_counter,:)))
    %check if any of the protein entries matahc the tails data
    [Position]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Protein_IDs_mito{unique_counter,terms_in_protein_group}, NTAILS.uniprot_accession_num, 'exact'));

    %test if
    if ~isempty(Position) & Checked_row_counter==0
        Overlap_array2_mito(Overlap_counter2_mito,1)=Protein_IDs_mito(unique_counter);
        Overlap_array2_mito(Overlap_counter2_mito,2)=NTAILS.uniprot_accession_num(Position(1));
        Overlap_counter2_mito=1+Overlap_counter2_mito;
        Checked_row_counter=1;
    end
     end
end

%determine the overlap between the mito and the cyto
Overlap_counter2_combined=1;

for unique_counter=2:length(Protein_IDs_mito)
    Checked_row_counter=0;    
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs_mito(unique_counter,:)))
    %check if any of the protein entries matahc the tails data
    [Position]=ind2sub(size(Protein_IDs), strmatch(Protein_IDs_mito{unique_counter,terms_in_protein_group}, Protein_IDs, 'exact'));

    %test if
    if ~isempty(Position) & Checked_row_counter==0
        Overlap_counter2_combined=1+Overlap_counter2_combined;
        Checked_row_counter=1;
    end
     end
end

%determine the overlap between the mito, the cyto and N-termini
Overlap_counter3_combined=1;

for unique_counter=2:length(Protein_IDs_mito)
    Checked_row_counter=0;    
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs_mito(unique_counter,:)))
    %check if any of the protein entries matahc the tails data
    [PositionA]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Protein_IDs_mito{unique_counter,terms_in_protein_group}, NTAILS.uniprot_accession_num, 'exact'));
    %check if any of the protein entries matahc the tails data
    [PositionB]=ind2sub(size(Protein_IDs), strmatch(Protein_IDs_mito{unique_counter,terms_in_protein_group}, Protein_IDs, 'exact'));

    %test if
    if ~isempty(PositionA) & ~isempty(PositionB) & Checked_row_counter==0
        Overlap_counter3_combined=1+Overlap_counter3_combined;
        Checked_row_counter=1;
    end
     end
end


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

%determine the Number of unique uniprot terms which we observe gaussians for and look up in the N-tails
%dataset
Unique_proteins_cyto=unique(Master_Gaussian_cyto.Proteinname);
Overlap_array_No_NaN1=cell(1,2);
Overlap_counter1_No_NaN=1;
Not_overlap_counter_No_Na=0;

for unique_counter=1:length(Unique_proteins_cyto)
    
    %look protein in major proein group file
    [row col]=ind2sub(size(Protein_IDs), strmatch(Unique_proteins_cyto{unique_counter}, Protein_IDs,'exact')); 
    
    Checked_row_counter=0;     
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs(row,:)))
        
    %check if any of the protein entries matahc the tails data
    [Position]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Protein_IDs{row,terms_in_protein_group}, NTAILS.uniprot_accession_num,'exact'));

    %test if
    if ~isempty(Position) & (~isnan(str2double(NTAILS.AverageNormalizedHvsL(Position(1)))) | ~isnan(str2double(NTAILS.AverageNormalizedMvsL(Position(1))))) & Checked_row_counter==0
        Overlap_array_No_NaN1(Overlap_counter1_No_NaN,1)=Unique_proteins_cyto(unique_counter);
        Overlap_array_No_NaN1(Overlap_counter1_No_NaN,2)=NTAILS.uniprot_accession_num(Position(1));
        Overlap_counter1_No_NaN=1+Overlap_counter1_No_NaN;
        Checked_row_counter=1;
    else
        Not_overlap_counter_No_Na=1+Not_overlap_counter_No_Na;
    end
    end
    
end

%look up in mito
Unique_proteins_mito=unique(Master_Gaussian_mito.Proteinname);
Overlap_array_No_NaN1_mito=cell(1,2);
Overlap_counter1_No_NaN_mito=1;
Not_overlap_counter_No_Na_mito=0;

for unique_counter=1:length(Unique_proteins_mito)
    
    %look protein in major proein group file
    [row col]=ind2sub(size(Protein_IDs_mito), strmatch(Unique_proteins_mito{unique_counter}, Protein_IDs_mito,'exact')); 
    
    Checked_row_counter=0;     
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs_mito(row,:)))
        
    %check if any of the protein entries matahc the tails data
    [Position]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Protein_IDs_mito{row,terms_in_protein_group}, NTAILS.uniprot_accession_num,'exact'));

    %test if
    if ~isempty(Position) & (~isnan(str2double(NTAILS.AverageNormalizedHvsL(Position(1)))) | ~isnan(str2double(NTAILS.AverageNormalizedMvsL(Position(1))))) & Checked_row_counter==0
        Overlap_array_No_NaN1_mito(Overlap_counter1_No_NaN_mito,1)=Unique_proteins_mito(unique_counter);
        Overlap_array_No_NaN1_mito(Overlap_counter1_No_NaN_mito,2)=NTAILS.uniprot_accession_num(Position(1));
        Overlap_counter1_No_NaN_mito=1+Overlap_counter1_No_NaN_mito;
        Checked_row_counter=1;
    else
        Not_overlap_counter_No_Na_mito=1+Not_overlap_counter_No_Na_mito;
    end
    end
    
end


%look up in mito compared to cyto
Overlap_counter1_No_NaN_mito2=1;
Not_overlap_counter_No_Na_mito2=0;

for unique_counter=1:length(Unique_proteins_mito)
    
    %look protein in major proein group file
    [rowA col]=ind2sub(size(Protein_IDs_mito), strmatch(Unique_proteins_mito{unique_counter}, Protein_IDs_mito,'exact')); 
    
    Checked_row_counter=0;     
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs_mito(rowA,:)))
        
    %check if any of the protein entries matahc the tails data
    [rowB col]=ind2sub(size(Protein_IDs), strmatch(Protein_IDs_mito{rowA,terms_in_protein_group}, Protein_IDs,'exact'));

    %test if
    if ~isempty(rowB) & Checked_row_counter==0
    %check if any of the protein entries matahc the tails data
    [Position]=ind2sub(size(Unique_proteins_cyto), strmatch(Protein_IDs{rowB,1}, Unique_proteins_cyto,'exact'));
    if ~isempty(Position) & Checked_row_counter==0
        Overlap_counter1_No_NaN_mito2=1+Overlap_counter1_No_NaN_mito2;
        Checked_row_counter=1;
    else
        Not_overlap_counter_No_Na_mito=1+Not_overlap_counter_No_Na_mito;
    end
    end
    end
    
end

%look up in mito compared to cyto
Overlap_counter1_Quantfied_all=1;

for unique_counter=1:length(Unique_proteins_mito)
    
    %look protein in major proein group file
    [rowA col]=ind2sub(size(Protein_IDs_mito), strmatch(Unique_proteins_mito{unique_counter}, Protein_IDs_mito,'exact')); 
    
    Checked_row_counter=0;     
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs_mito(rowA,:)))
        
    %check if any of the protein entries matahc the tails data
    [rowB col]=ind2sub(size(Protein_IDs), strmatch(Protein_IDs_mito{rowA,terms_in_protein_group}, Protein_IDs,'exact'));
    
    %test if
    if ~isempty(rowB) & Checked_row_counter==0
    
    %check if any of the protein entries matahc the tails data
    [PositionA]=ind2sub(size(Unique_proteins_cyto), strmatch(Protein_IDs{rowB,1}, Unique_proteins_cyto,'exact'));
    
    %check if any of the protein entries matahc the tails data
    [PositionB]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Unique_proteins_cyto{PositionA}, NTAILS.uniprot_accession_num,'exact'));

    
    if ~isempty(PositionB)& (~isnan(str2double(NTAILS.AverageNormalizedHvsL(PositionB(1)))) & ~isnan(str2double(NTAILS.AverageNormalizedMvsL(PositionB(1))))) & Checked_row_counter==0
        Overlap_counter1_Quantfied_all=1+Overlap_counter1_Quantfied_all;
        Checked_row_counter=1;
    end
    end
    end
    
end


%determine the Number of unique uniprot terms which we observe gaussians for and look up in the N-tails
%dataset
Overlap_array_No_NaN2=cell(length(NTAILS.uniprot_accession_num),2);
Overlap_counter_No_NaN2=1;

for unique_counter=2:length(Protein_IDs)
    Rowed_checked=0;    
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs(unique_counter,:)))
    %check if any of the protein entries matahc the tails data
    [Position]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Protein_IDs{unique_counter,terms_in_protein_group}, NTAILS.uniprot_accession_num,'exact'));

    %test if
    if ~isempty(Position) & (~isnan(str2double(NTAILS.AverageNormalizedHvsL(Position(1)))) | ~isnan(str2double(NTAILS.AverageNormalizedMvsL(Position(1))))) & Rowed_checked==0
        Overlap_array_No_NaN2(Position(1),1)=Protein_IDs(unique_counter);
        Overlap_array_No_NaN2(Position(1),2)=NTAILS.uniprot_accession_num(Position(1));
        Overlap_counter_No_NaN2=1+Overlap_counter_No_NaN2;
        Rowed_checked=1;
    end
     end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%determine N tail proteins with values
Number_unique_Ntails_proteins=NaN(length(Unique_N_tails_uniprot),2);

for unique_tail_counter=1:length(Unique_N_tails_uniprot)

    [Position1]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Unique_N_tails_uniprot{unique_tail_counter}, NTAILS.uniprot_accession_num,'exact'));
    
    TempA=zeros(1,1);
    TempB=zeros(1,1);
    for temp_counter=1:length(Position1)
    TempA(temp_counter)=str2double(NTAILS.AverageNormalizedMvsL{Position1(temp_counter)});
    TempB(temp_counter)=str2double(NTAILS.AverageNormalizedHvsL{Position1(temp_counter)});
    end
    
    Number_unique_Ntails_proteins(unique_tail_counter,1)=nanmean(TempA);
    Number_unique_Ntails_proteins(unique_tail_counter,2)=nanmean(TempB);
    
end

%Remove row where both values are NaN
[rows_A] = find(isnan(Number_unique_Ntails_proteins(:,1)));
[rows_B] = find(isnan(Number_unique_Ntails_proteins(:,1)));
merged_row=intersect(rows_A,rows_B);
Number_unique_Ntails_proteins(merged_row,:) = [];

%Determine if there is a correlation with proteins that change between the
%two type of MS data
%work out unique proteins in Diversity array
Unique_proteins_Diversity=unique(diversity_in_peptide.Protein);

Overlap_array_for_correlation=zeros(1,3);
Proteins_which_overlap=cell(1,2);
Number_of_terminin_per_proteins=zeros(1,1);
Signficant_termini=cell(1,2);
counter=1;


for Protein_counter1=1:length(Unique_proteins_Diversity)
    
    %Find protein in array 
    %check if any of the protein entries matahc the tails data
    [Position1]=ind2sub(size(diversity_in_peptide.Protein), strmatch(Unique_proteins_Diversity{Protein_counter1}, diversity_in_peptide.Protein,'exact'));
    Temp_value=nanmean(diversity_in_peptide.Evenness(Position1));
    
    %Find protein in list of protein groups
    [PositionA, PositionB]=ind2sub(size(Protein_IDs), strmatch(Unique_proteins_Diversity{Protein_counter1}, Protein_IDs,'exact'));
    
    Checked_row_counter=0;
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs(PositionA,:)))
    
    %check if any of the protein entries matahc the tails data
    PositionC=[];
    [PositionC]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Protein_IDs{PositionA,terms_in_protein_group}, NTAILS.uniprot_accession_num));

    %test if
    if ~isempty(PositionC) & Checked_row_counter==0
    
    %determine the unique peptide
    A_1=unique(NTAILS.Peptidesequence(PositionC));
    B_1=NTAILS.Peptidesequence(PositionC);
    %Find position of unique terms
    [~,unique_temp]=ismember(A_1,B_1);
    
    Temp_value2_HvsL=NaN(1,1);
    Temp_value3_MvsL=NaN(1,1);
    Signficant_terminiMWW=NaN(1,1);
    Signficant_termini_ttest=NaN(1,1);
    
    for unique_counter=1:length(unique_temp)
    
    %Determine average peptide measurements
    Temp_value2_HvsL(unique_counter)=str2double(NTAILS.AverageNormalizedHvsL(PositionC(unique_temp(unique_counter))));
    Temp_value3_MvsL(unique_counter)=str2double(NTAILS.AverageNormalizedMvsL(PositionC(unique_temp(unique_counter))));
    
    %check if any of the termini show signficant chnages
    if str2double(NTAILS.BenjaminiHochbergcorrectedpvalueMvsLMWW(PositionC(unique_temp(unique_counter))))<0.05 |...
       str2double(NTAILS.BenjaminiHochbergcorrectedpvalueHvsLMWW(PositionC(unique_temp(unique_counter))))<0.05 
       Signficant_terminiMWW=1;    
    end
    
    %check if any of the termini show signficant chnages
    if str2double(NTAILS.BenjaminiHochbergcorrectedpvalueMvsLttest(PositionC(unique_temp(unique_counter))))<0.05 |...
       str2double(NTAILS.BenjaminiHochbergcorrectedpvalueHvsLttest(PositionC(unique_temp(unique_counter))))<0.05 
       Signficant_termini_ttest=1;    
    end
    
    end
    
    %Write out values
    Proteins_which_overlap(counter,1)=Unique_proteins_Diversity(Protein_counter1);
    Proteins_which_overlap(counter,2)=NTAILS.uniprot_accession_num(PositionC(1));
    Overlap_array_for_correlation(counter,1)=Temp_value;
    Overlap_array_for_correlation(counter,3)=nanmean(Temp_value2_HvsL);
    Overlap_array_for_correlation(counter,2)=nanmean(Temp_value3_MvsL);
    Number_of_terminin_per_proteins(counter)=length(unique_temp);
    
    if Signficant_terminiMWW==1
    Signficant_termini(counter,1)={'+'};
    else
    Signficant_termini(counter,1)={''};
    end
    
    if Signficant_termini_ttest==1
    Signficant_termini(counter,2)={'+'};
    else
    Signficant_termini(counter,2)={''};
    end    
    
    counter=counter+1;
    Checked_row_counter=1;
        
    end
    end
end
    
%Combine Persaus Gaussian table and N-tails together
Output_Ntails_sign=cell(length(Perseus_enrichment_Protein_level.Proteinname),2);

for persaus_counter=2:length(Perseus_enrichment_Protein_level.Proteinname)

    %look protein in major proein group file
    [row col]=ind2sub(size(Protein_IDs), strmatch(Perseus_enrichment_Protein_level.Proteinname{persaus_counter}, Protein_IDs,'exact')); 
    
    Checked_row_counter=0;
    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs(row,:)))
        
    Position=[];   
    %check if any of the protein entries matahc the tails data
    [Position]=ind2sub(size(Proteins_which_overlap(:,2)), strmatch(Protein_IDs{row,terms_in_protein_group}, Proteins_which_overlap(:,2), 'exact'));

    %test if
    if ~isempty(Position) & Checked_row_counter==0

    Output_Ntails_sign(persaus_counter,1)=Signficant_termini(Position,1);
    Output_Ntails_sign(persaus_counter,2)=Signficant_termini(Position,2);
    Checked_row_counter=1;    
    end
    end  
    
end

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


%Write out table for enrichment analysis
Output_name=strcat('Persaus_Protein_level_Gaussian_and_Ntails.txt');

fid = fopen(Output_name,'w');
fprintf (fid,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t\n',...
'Protein name',...
'Increase AND Decreases (>= two replicates)',...
'Increase (>= two replicates)',...
'Decrease (>= two replicates)',...
'Gaussian fitted with change below adjusted_pvalue(bonferroni correction) (ttest)',...
'Gaussian fitted with change below adjusted_pvalue(bonferroni correction) (MWW)',...
'N-tail change below adjusted_pvalue(Ben Hochberg) (ttest)',....
'N-tail change below adjusted_pvalue(Ben Hochberg) (MWW)'); %Write Header
for counter_output=2:length(Perseus_enrichment_Protein_level.Proteinname)
fprintf(fid,'%s	',Perseus_enrichment_Protein_level.Proteinname{counter_output});
fprintf(fid,'%s	',Perseus_enrichment_Protein_level.IncreaseANDDecreasestworeplicates{counter_output});
fprintf(fid,'%s	',Perseus_enrichment_Protein_level.Increasetworeplicates{counter_output});
fprintf(fid,'%s	',Perseus_enrichment_Protein_level.Decreasetworeplicates{counter_output});
fprintf(fid,'%s	',Perseus_enrichment_Protein_level.ttest_pvaluebonferronicorrectio{counter_output});
fprintf(fid,'%s	',Perseus_enrichment_Protein_level.MWW_pvaluebonferronicorrecti1{counter_output});
fprintf(fid,'%s	',Output_Ntails_sign{counter_output,2});
fprintf(fid,'%s	',Output_Ntails_sign{counter_output,1});
fprintf(fid,'\n');
end
fclose(fid);

% Group based on individual Termini
Overlap_array_for_correlation_termini=zeros(1,3);
Proteins_which_overlap_termini=cell(1,1);
Ntermini_sequence=cell(1,1);
Number_of_terminin=zeros(1,1);
counter=1;

for Protein_counter1=1:length(Unique_proteins_Diversity)
    
    %Find protein in array 
    %check if any of the protein entries matahc the tails data
    [Position1]=ind2sub(size(diversity_in_peptide.Protein), strmatch(Unique_proteins_Diversity{Protein_counter1}, diversity_in_peptide.Protein,'exact'));
    Temp_value=nanmean(diversity_in_peptide.Evenness(Position1));
    
    %Find protein in list of protein groups
    [PositionA, PositionB]=ind2sub(size(Protein_IDs), strmatch(Unique_proteins_Diversity{Protein_counter1}, Protein_IDs,'exact'));

    for terms_in_protein_group=1:sum(~cellfun('isempty',Protein_IDs(PositionA,:)))
    
    %reset value
    Position2=[];
    %check if any of the protein entries matahc the tails data
    [Position2]=ind2sub(size(NTAILS.uniprot_accession_num), strmatch(Protein_IDs{PositionA,terms_in_protein_group}, NTAILS.uniprot_accession_num));

    %test if
    if ~isempty(Position2)
        
    %determine the unique peptide
    A_1=unique(NTAILS.Peptidesequence(Position2));
    B_1=NTAILS.Peptidesequence(Position2);
    %Find position of unique terms
    [~,unique_temp]=ismember(A_1,B_1);
    
    
    %Determine average peptide measurements
    for converting_counter=1:length(unique_temp)
    Temp_value2_HvsL=str2double(NTAILS.AverageNormalizedHvsL(Position2(unique_temp(converting_counter))));
    Temp_value3_MvsL=str2double(NTAILS.AverageNormalizedMvsL(Position2(unique_temp(converting_counter))));


    %Write out values
    Proteins_which_overlap_termini(counter)=Unique_proteins_Diversity(Protein_counter1);
    Ntermini_sequence(counter)=NTAILS.PeptideSequencewithmodifications(Position2(converting_counter));
    Overlap_array_for_correlation_termini(counter,1)=Temp_value;
    Overlap_array_for_correlation_termini(counter,3)=Temp_value2_HvsL;
    Overlap_array_for_correlation_termini(counter,2)=Temp_value3_MvsL;
    counter=counter+1;
    end
    %save number ot termini pre proteins
    Number_of_terminin(counter)=length(unique_temp);
    end
    end

end

%Calicualte correlation on the protein level
%Create list without NaN
Overlap_array_for_correlation_NoNAN=Overlap_array_for_correlation;

[rows, cols] = find(isnan(Overlap_array_for_correlation_NoNAN)); 
Overlap_array_for_correlation_NoNAN(unique(rows),:) = [];

%Determine correlation
[RHO_MvsL,PVAL_MvsL] = corr(Overlap_array_for_correlation_NoNAN(:,1),Overlap_array_for_correlation_NoNAN(:,2),'Type','Spearman');
[RHO_HvsL,PVAL_HvsL] = corr(Overlap_array_for_correlation_NoNAN(:,1),Overlap_array_for_correlation_NoNAN(:,3),'Type','Spearman');


%Calicualte correlation on the protein level
%Create list without NaN
Overlap_array_for_correlation_termini_NoNAN=Overlap_array_for_correlation_termini;

[rows, cols] = find(isnan(Overlap_array_for_correlation_termini_NoNAN)); 
Overlap_array_for_correlation_termini_NoNAN(unique(rows),:) = [];

%Determine correlation
[RHO_termini_MvsL,PVAL_termini_MvsL] = corr(Overlap_array_for_correlation_termini_NoNAN(:,1),Overlap_array_for_correlation_termini_NoNAN(:,2),'Type','Spearman');
[RHO_termini_HvsL,PVAL_termini_HvsL] = corr(Overlap_array_for_correlation_termini_NoNAN(:,1),Overlap_array_for_correlation_termini_NoNAN(:,3),'Type','Spearman');

F2B=figure
scatter(log2(Overlap_array_for_correlation_NoNAN(:,1)),log2(Overlap_array_for_correlation_NoNAN(:,2)),'filled')
hold on
scatter(log2(Overlap_array_for_correlation_NoNAN(:,1)),log2(Overlap_array_for_correlation_NoNAN(:,3)),'filled')
xlabel('log2 of Unevenness within peptide measurement');
ylabel('log2 Average N-TAILS ratio compared to untreated cells');
title('Peptide uneveness versus average N-TAILS ratio (Averaged per protein)');
legend('4 hour Fas treatment', '4 hour Fas treatment with inhibitor','FontSize',6, 'Location', 'southoutside')
Output_name8=strcat('Scatter_plot_uneveness_versus_ratio_proteins.pdf');
print(F2B,'-dpdf','-r600', Output_name8);


F4=figure
scatter(log2(Overlap_array_for_correlation(:,1)),log2(Overlap_array_for_correlation(:,2)),'filled')
hold on
scatter(log2(Overlap_array_for_correlation(:,1)),log2(Overlap_array_for_correlation(:,3)),'filled')
xlabel('log2 of Unevenness within peptide measurement');
ylabel('log2 Average N-TAILS ratio compared to untreated cells');
title('Peptide uneveness versus average N-TAILS ratio (Averaged per protein)');
legend('4 hour Fas treatment', '4 hour Fas treatment with inhibitor','FontSize',6, 'Location', 'southoutside')
Output_name8=strcat('Scatter_plot_uneveness_versus_ratio.pdf');
print(F4,'-dpdf','-r600', Output_name8);


F5=figure
scatter(log2(Overlap_array_for_correlation_termini(:,1)),log2(Overlap_array_for_correlation_termini(:,2)),'filled')
hold on
scatter(log2(Overlap_array_for_correlation_termini(:,1)),log2(Overlap_array_for_correlation_termini(:,3)),'filled')
xlabel('log2 of Unevenness within peptide measurement');
ylabel('Log2 Average N-TAILS ratio compared to untreated cells');
title('Peptide uneveness versus average N-TAILS ratio (per termini)');
legend('4 hour Fas treatment', '4 hour Fas treatment with inhibitor','FontSize',6, 'Location', 'southoutside')
Output_name8=strcat('Scatter_plot_uneveness_versus_ratio.pdf');
print(F5,'-dpdf','-r600', Output_name8);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Determine the average diversity per termini
max_number_termini=max(Number_of_terminin_per_proteins);
Average_eveness_observed=zeros(3,max(Number_of_terminin_per_proteins));

for max_counter=1:max_number_termini
    %Find positions of termini equal to value being assessed
    [sort_orderA sort_orderB]=find(Number_of_terminin_per_proteins==max_counter);
    
    if length(sort_orderB)>5
    try
    Average_eveness_observed(1,max_counter)=nanmean(Overlap_array_for_correlation(sort_orderB,1));
    Average_eveness_observed(2,max_counter)=max(Overlap_array_for_correlation(sort_orderB,1));
    Average_eveness_observed(3,max_counter)=nanstd(Overlap_array_for_correlation(sort_orderB,1));
    catch
    Average_eveness_observed(1,max_counter)=0;
    Average_eveness_observed(2,max_counter)=0;  
    Average_eveness_observed(3,max_counter)=0; 
    end
    else 
    Average_eveness_observed(1,max_counter)=NaN;
    Average_eveness_observed(2,max_counter)=NaN;  
    Average_eveness_observed(3,max_counter)=NaN; 
    end
    
end

total_mean_evenness=nanmean(Overlap_array_for_correlation(:,1));

F5=figure
scatter(Number_of_terminin_per_proteins,log10(Overlap_array_for_correlation(:,1)),'filled')
hold on
scatter(1:max(Number_of_terminin_per_proteins),log10(Average_eveness_observed(1,:)),'filled')
ylim([0 log10(max(Average_eveness_observed(2,8)))])
xlim([0 8.5])
ylabel('Average Unevenness within proteins measurement');
xlabel('Number of N-termini detected');
title('Protein uneveness versus average N-terminal peptide detected');
hold on
plot([1 10],[log10(total_mean_evenness) log10(total_mean_evenness)],':k');
legend('mean evenness measurements', 'Standard deviation of evenness','mean uneveness','FontSize',6, 'Location', 'eastoutside')
Output_name8=strcat('Scatter_plot_uneveness_versus_N_termini_number_all_data_point.pdf');
print(F5,'-dpdf','-r600', Output_name8);


F5=figure
scatter(1:max(Number_of_terminin_per_proteins),Average_eveness_observed(1,:),'filled')
hold on
plot([1 10],[total_mean_evenness total_mean_evenness],':k');
ylabel('Average Unevenness within proteins measurement');
xlabel('Number of N-termini detected');
title('Protein uneveness versus average N-terminal peptide detected');
legend('mean evenness measurements', 'Standard deviation of evenness','mean uneveness','FontSize',6, 'Location', 'eastoutside')
Output_name8=strcat('Scatter_plot_uneveness_versus_N_termini_number_all_data_point.pdf');
print(F5,'-dpdf','-r600', Output_name8);


Output_name1=strcat('NTAILS_and peptide_unevenness_analysis_protein_level.txt');

fid = fopen(Output_name1,'w');
fprintf (fid,'%s\t%s\t%s\t%s\t\n',...
'Protein observed in both NTAILS and interactome',...
'Average unevenness observed within interactome',...
'N-TAILS ratio observed in 4 hour Fas treatment',...
'N-TAILS ratio observed in 4 hour Fas treatment with inhibitor'); %Write Header
for counter_output=1:length(Proteins_which_overlap)
fprintf(fid,'%s	',Proteins_which_overlap{counter_output});
fprintf(fid,'%f	',Overlap_array_for_correlation(counter_output,1));
fprintf(fid,'%f	',Overlap_array_for_correlation(counter_output,2));
fprintf(fid,'%f	',Overlap_array_for_correlation(counter_output,3));
fprintf(fid,'\n');
end
fclose(fid);


Output_name2=strcat('Summary_NTAILS_and peptide_unevenness_analysis_Protein_level.txt');
[h,p1,ci,stats] = ttest2(Overlap_array_for_correlation(:,2),Overlap_array_for_correlation(:,3));

fid = fopen(Output_name2,'w');
fprintf (fid,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t',...
'Number of termini observed for protein identified within the interactome',...
'Mean Uneveness of protein with N-termini observed',...
'Mean N-TAILS ratio observed in 4 hour Fas treatment',...
'Mean N-TAILS ratio observed in 4 hour Fas treatment with inhibitor',...
'Median N-TAILS ratio observed in 4 hour Fas treatment',...
'Median N-TAILS ratio observed in 4 hour Fas treatment with inhibitor',...
'two sideded t-test p-value','Spearman correlation Rho (MvsL)',...
'Spearman correlation p-value (MvsL)''Spearman correlation Rho (HvsL)',...
'Spearman correlation p-value (HvsL)'); %Write Header
fprintf(fid,'\n');
fprintf(fid,'%f	',length(Overlap_array_for_correlation));
fprintf(fid,'%f	',nanmean(Overlap_array_for_correlation(:,1)));
fprintf(fid,'%f	',nanmean(Overlap_array_for_correlation(:,2)));
fprintf(fid,'%f	',nanmean(Overlap_array_for_correlation(:,3)));
fprintf(fid,'%f	',nanmedian(Overlap_array_for_correlation(:,2)));
fprintf(fid,'%f	',nanmedian(Overlap_array_for_correlation(:,3)));
fprintf(fid,'%f	',p1);
fprintf(fid,'%f	',RHO_MvsL);
fprintf(fid,'%f	',PVAL_MvsL);
fprintf(fid,'%f	',RHO_HvsL);
fprintf(fid,'%f	',PVAL_HvsL);
fclose(fid);


Output_name3=strcat('NTAILS_and peptide_unevenness_analysis_Ntermini_level.txt');

fid = fopen(Output_name3,'w');
fprintf (fid,'%s\t%s\t%s\t%s\t%s\t\n',...
'Protein observed in both NTAILS and interactome',...
'N_termini observed',...
'Average unevenness observed within interactome',...
'N-TAILS ratio observed in 4 hour Fas treatment',...
'N-TAILS ratio observed in 4 hour Fas treatment with inhibitor'); %Write Header
for counter_output=1:length(Proteins_which_overlap_termini)
fprintf(fid,'%s	',Proteins_which_overlap_termini{counter_output});
fprintf(fid,'%s	',Ntermini_sequence{counter_output});
fprintf(fid,'%f	',Overlap_array_for_correlation_termini(counter_output,1));
fprintf(fid,'%f	',Overlap_array_for_correlation_termini(counter_output,2));
fprintf(fid,'%f	',Overlap_array_for_correlation_termini(counter_output,3));
fprintf(fid,'\n');
end
fclose(fid);

Output_name4=strcat('Summary_NTAILS_and peptide_unevenness_analysis_Ntermini_level.txt');
[h,p,ci,stats] = ttest2(Overlap_array_for_correlation_termini(:,2),Overlap_array_for_correlation_termini(:,3));

fid = fopen(Output_name4,'w');
fprintf (fid,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t',...
'Number of termini observed for protein identified within the interactome',...
'Mean Uneveness of protein with N-termini observed',...
'Mean N-TAILS ratio observed in 4 hour Fas treatment',...
'Mean N-TAILS ratio observed in 4 hour Fas treatment with inhibitor',...
'Median N-TAILS ratio observed in 4 hour Fas treatment',...
'Median N-TAILS ratio observed in 4 hour Fas treatment with inhibitor',...
'two sideded t-test p-value','Spearman correlation Rho (MvsL)',...
'Spearman correlation p-value (MvsL)''Spearman correlation Rho (HvsL)',...
'Spearman correlation p-value (HvsL)'); %Write Header
fprintf(fid,'\n');
fprintf(fid,'%f	',length(Overlap_array_for_correlation_termini));
fprintf(fid,'%f	',nanmean(Overlap_array_for_correlation_termini(:,1)));
fprintf(fid,'%f	',nanmean(Overlap_array_for_correlation_termini(:,2)));
fprintf(fid,'%f	',nanmean(Overlap_array_for_correlation_termini(:,3)));
fprintf(fid,'%f	',nanmedian(Overlap_array_for_correlation_termini(:,2)));
fprintf(fid,'%f	',nanmedian(Overlap_array_for_correlation_termini(:,3)));
fprintf(fid,'%f	',p);
fprintf(fid,'%f	',RHO_termini_MvsL);
fprintf(fid,'%f	',PVAL_termini_MvsL);
fprintf(fid,'%f	',RHO_termini_HvsL);
fprintf(fid,'%f	',PVAL_termini_HvsL);
fclose(fid);


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

two_D_array=zeros(2,1);
two_D_arrayB=zeros(2,1);
two_D_arrayC=zeros(2,1);
two_D_array_justAsp=zeros(2,1);
two_D_array_justArg=zeros(2,1);
Gaus_counterA=1;
Gaus_counterB=1;
Gaus_counterC=1;
Gaus_counterD=1;
Gaus_counterE=1;
temp_counter=1;

for max_counter=1:length(Proteins_which_overlap(:,1))

%Find how many termini were detected
[Gaussian_list_position]=ind2sub(size(Imported_seperated_cyto(:,1)),...
strmatch(Proteins_which_overlap{max_counter,1}, Imported_seperated_cyto(:,1), 'exact')); 

%Find how many termini were detected
[Position_in_NTAILS_list]=ind2sub(size(NTAILS.uniprot_accession_num(:,1)),...
strmatch(Proteins_which_overlap{max_counter,2}, NTAILS.uniprot_accession_num(:,1), 'exact'));

for Gaus_number_counte=1:length(Gaussian_list_position)
for Ntail_number_counte=1:length(Position_in_NTAILS_list)
    if ~isequal(NTAILS.Previousaminoacid(Position_in_NTAILS_list(Ntail_number_counte)),{'D'})
    two_D_array(1,Gaus_counterA)=str2num(Imported_seperated_cyto{Gaussian_list_position(Gaus_number_counte),18});
    two_D_array(2,Gaus_counterA)=log2(str2num(NTAILS.AverageNormalizedHvsL{Position_in_NTAILS_list(Ntail_number_counte)}));
    Gaus_counterA=Gaus_counterA+1;
    else
    two_D_array_justAsp(1,Gaus_counterB)=str2num(Imported_seperated_cyto{Gaussian_list_position(Gaus_number_counte),18});
    two_D_array_justAsp(2,Gaus_counterB)=log2(str2num(NTAILS.AverageNormalizedHvsL{Position_in_NTAILS_list(Ntail_number_counte)}));
    Gaus_counterB=Gaus_counterB+1;
    end
end
end

for Gaus_number_counte=1:length(Gaussian_list_position)
for Ntail_number_counte=1:length(Position_in_NTAILS_list)
    if isequal(NTAILS.Previousaminoacid(Position_in_NTAILS_list(Ntail_number_counte)),{'R'})
    two_D_array_justArg(1,Gaus_counterE)=str2num(Imported_seperated_cyto{Gaussian_list_position(Gaus_number_counte),18});
    two_D_array_justArg(2,Gaus_counterE)=log2(str2num(NTAILS.AverageNormalizedHvsL{Position_in_NTAILS_list(Ntail_number_counte)}));
    Gaus_counterE=Gaus_counterE+1;
    end
end
end

TempA=zeros(2,1);
TempB=zeros(2,2);

%Only record Max and Min values
for Gaus_number_counte=1:length(Gaussian_list_position)
for Ntail_number_counte=1:length(Position_in_NTAILS_list)

    TempA(1,temp_counter)=str2num(Imported_seperated_cyto{Gaussian_list_position(Gaus_number_counte),18});
    TempA(2,temp_counter)=log2(str2num(NTAILS.AverageNormalizedHvsL{Position_in_NTAILS_list(Ntail_number_counte)}));
    temp_counter=temp_counter+1;

end
end
TempA(TempA==0)=NaN;

two_D_arrayB(1,Gaus_counterC)=max(TempA(1,:));
two_D_arrayB(2,Gaus_counterC)=max(TempA(2,:));
Gaus_counterC=Gaus_counterC+1;

two_D_arrayC(1,Gaus_counterD)=min(TempA(1,:));
two_D_arrayC(2,Gaus_counterD)=min(TempA(2,:));
Gaus_counterD=Gaus_counterD+1;

end

two_D_array(two_D_array==0)=NaN;
two_D_array_justAsp(two_D_array_justAsp==0)=NaN;
two_D_array_justArg(two_D_array_justArg==0)=NaN;
two_D_arrayB(two_D_arrayB==0)=NaN;
two_D_arrayC(two_D_arrayC==0)=NaN;

%Determine mean and median, need to plot these also add lines at zero
Mean_N_termini_arrayB=nanmean(two_D_arrayB(2,:));
Median_N_termini_arrayB=nanmedian(two_D_arrayB(2,:));
Mean_N_termini_arrayC=nanmean(two_D_arrayC(2,:));
Median_N_termini_arrayC=nanmedian(two_D_arrayC(2,:));
Mean_N_termini=nanmean(two_D_array(2,:));
Median_N_termini=nanmedian(two_D_array(2,:));
Mean_Guas_arrayB=nanmean(two_D_arrayB(1,:));
Median_Guas_arrayB=nanmedian(two_D_arrayB(1,:));
Mean_Guas_arrayC=nanmean(two_D_arrayC(1,:));
Median_Guas_arrayC=nanmedian(two_D_arrayC(1,:));
Mean_Guas=nanmean(two_D_array(1,:));
Median_Guas=nanmedian(two_D_array(1,:));

%Median of just asp cleavages
Median_Asp_p1=nanmedian(two_D_array_justAsp(2,:));
Median_Asp_Gaus=nanmedian(two_D_array_justAsp(1,:));
%Median of Arg p1 peptide
Median_Arg_p1=nanmedian(two_D_array_justArg(2,:));
Median_Arg_Gaus=nanmedian(two_D_array_justArg(1,:));

%Create text showing median
Text1A=strcat('Median Gaussian fold change : ', mat2str(round(Median_Guas*100)/100));
Text2A=strcat('Median N-termini ratio', mat2str(round(Median_N_termini*100)/100));
TextAsp_C=strcat('Median Asp P1 N-termini ratio :', mat2str(round(Median_Asp_p1*100)/100));
TextArg_C=strcat('Median Arg P1 N-termini ratio :', mat2str(round(Median_Arg_p1*100)/100));
Text1B=strcat('Median Gaussian fold change : ', mat2str(round(Median_Guas_arrayB*100)/100));
Text2B=strcat('Median N-termini ratio', mat2str(round(Median_N_termini_arrayB*100)/100));
Text1D=strcat('Median Gaussian fold change : ', mat2str(round(Median_Guas_arrayC*100)/100));
Text2D=strcat('Median N-termini ratio', mat2str(round(Median_N_termini_arrayC*100)/100));

%Create scatter plot
scatter_plot_1=figure
scatter(two_D_array(1,:),two_D_array(2,:),'filled')
hold on
scatter(two_D_array_justAsp(1,:),two_D_array_justAsp(2,:),'filled')
%hold on
%scatter(two_D_array_justArg(1,:),two_D_array_justArg(2,:),'filled')
ylabel('N-Tail log2 ratios');
xlabel('Gaussian log2 ratios');
%hold on
%plot([-2:0.2:5],repmat(Median_N_termini,length([-2:0.2:5])),'-.r','linewidth', 1);
%hold on
%plot(repmat(Median_Guas,length([-4:0.2:14])),[-4:0.2:14],'-.b','linewidth', 1);
hold on
plot(repmat(0,length([-4:0.2:14])),[-4:0.2:14],'-.k','linewidth', 1);
hold on
plot([-2:0.2:5],repmat(0,length([-2:0.2:5])),'-.k','linewidth', 1);
legend('All P1 position','Aspartic acid P1 position',sprintf(Text1A),sprintf(Text2A),sprintf(TextAsp_C),sprintf(TextArg_C));
title('Comparsion of All N-tail ratio to all gaussian fold change combinations');
print(scatter_plot_1,'-dpdf','-r600', 'Scatterplot_log2_Ntails_vs_Gaus_changes.pdf');

scatter_plot_1=figure
scatter(two_D_array(1,:),two_D_array(2,:),'filled')
ylabel('N-Tail log2 ratios');
xlabel('Gaussian log2 ratios');
hold on
plot(repmat(0,length([-4:0.2:14])),[-4:0.2:14],'-.k','linewidth', 1);
hold on
plot([-2:0.2:5],repmat(0,length([-2:0.2:5])),'-.k','linewidth', 1);
title('Comparsion of All N-tail ratio to all gaussian fold change combinations');
print(scatter_plot_1,'-dpdf','-r600', 'Scatterplot_log2_Ntails_vs_Gaus_changes_version2.pdf');


%Create scatter plot
scatter_plot_1=figure
%scatter(two_D_arrayC(1,:),two_D_arrayC(2,:),'filled')
%hold on
scatter(two_D_arrayB(1,:),two_D_arrayB(2,:),'filled')
ylabel('N-Tail log2 ratios');
xlabel('Gaussian log2 ratios');
hold on
plot([-2:0.2:5],repmat(Median_N_termini_arrayB,length([-2:0.2:5])),'-.r','linewidth', 1);
hold on
plot(repmat(Median_Guas_arrayB,length([-4:0.2:14])),[-4:0.2:14],'-.b','linewidth', 1);
legend('All P1 position',sprintf(Text1B),sprintf(Text2B));
title('Comparsion of maximum N-tail ratio to maximum gaussian fold change');

print(scatter_plot_1,'-dpdf','-r600', 'Scatterplot_log2_Ntails_vs_Gaus_changes_max_combinations.pdf');


scatter_plot_1=figure
%scatter(two_D_arrayC(1,:),two_D_arrayC(2,:),'filled')
%hold on
scatter(two_D_arrayC(1,:),two_D_arrayB(2,:),'filled')
ylabel('N-Tail log2 ratios');
xlabel('Gaussian log2 ratios');
hold on
plot([-2:0.2:2.5],repmat(Median_N_termini_arrayC,length([-2:0.2:2.5])),'-.r','linewidth', 1);
hold on
plot(repmat(Median_Guas_arrayC,length([-4:0.2:14])),[-4:0.2:14],'-.b','linewidth', 1);
legend('All P1 position',sprintf(Text1D),sprintf(Text2D));
title('Comparsion of maximum N-tail ratio to minimum gaussian fold change');

print(scatter_plot_1,'-dpdf','-r600', 'Scatterplot_log2_Ntails_vs_Gaus_changes_max_combinations.pdf');

%remove nan
temp_two_D_arrayC=two_D_arrayC;
temp_two_D_arrayB=two_D_arrayB;

[col, rowA]=find(~isnan(temp_two_D_arrayB(1,:))==1);
[col, rowB]=find(~isnan(temp_two_D_arrayB(2,:))==1);
[col, rowC]=find(~isnan(temp_two_D_arrayC(1,:))==1);
[col, rowD]=find(~isnan(temp_two_D_arrayC(2,:))==1);

tempA=intersect(rowB,rowA);
tempB=intersect(rowC,rowD);
tempC=intersect(tempB,tempA);

[RHO_A,PVAL_A] = corr(two_D_arrayC(1,tempC)',two_D_arrayB(2,tempC)','Type','Spearman');
[RHO_B,PVAL_B] = corr(two_D_arrayC(1,tempC)',two_D_arrayB(2,tempC)','Type','Pearson');

%Write out table of distrubtion of N-termini peptide
fid = fopen('Correlation_between_N-termini_and_Gaus_observed_to_change.txt','w');
fprintf (fid,'%s\t%s\t%s\t',...
'Correlation test',...
'Correlation value',...
'P-value'); %Write Header
fprintf(fid,'\n');
fprintf(fid,'%s	','Spearman');    
fprintf(fid,'%f	',RHO_A(1));
fprintf(fid,'%f	',PVAL_A(1));
fprintf(fid,'\n');
fprintf(fid,'%s	','Pearson');    
fprintf(fid,'%f	',RHO_B(1));
fprintf(fid,'%f	',PVAL_B(1));
fprintf(fid,'\n');
fclose(fid);

%%% correlation to all combinations

%remove nan
temp_two_D_arrayC=two_D_array;
temp_two_D_arrayB=two_D_array;

[col, rowA]=find(~isnan(temp_two_D_arrayB(1,:))==1);
[col, rowB]=find(~isnan(temp_two_D_arrayB(2,:))==1);
[col, rowC]=find(~isnan(temp_two_D_arrayC(1,:))==1);
[col, rowD]=find(~isnan(temp_two_D_arrayC(2,:))==1);

tempA=intersect(rowB,rowA);
tempB=intersect(rowC,rowD);
tempC=intersect(tempB,tempA);

[RHO_A,PVAL_A] = corr(two_D_array(1,tempC)',two_D_array(2,tempC)','Type','Spearman');
[RHO_B,PVAL_B] = corr(two_D_array(1,tempC)',two_D_array(2,tempC)','Type','Pearson');

%Write out table of distrubtion of N-termini peptide
fid = fopen('Correlation_between_All_N-termini_and_Gaus_combinations_observed_to_change.txt','w');
fprintf (fid,'%s\t%s\t%s\t',...
'Correlation test',...
'Correlation value',...
'P-value'); %Write Header
fprintf(fid,'\n');
fprintf(fid,'%s	','Spearman');    
fprintf(fid,'%f	',RHO_A(1));
fprintf(fid,'%f	',PVAL_A(1));
fprintf(fid,'\n');
fprintf(fid,'%s	','Pearson');    
fprintf(fid,'%f	',RHO_B(1));
fprintf(fid,'%f	',PVAL_B(1));
fprintf(fid,'\n');
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Determine number of interacomte changes which correlate with N-terminal

two_D_array2=two_D_array(:,sum(~isnan(two_D_array))==2);
Counter_chnages=0;

for counter_1=1:length(two_D_array2)
    if (two_D_array2(1,counter_1)>1 | two_D_array2(1,counter_1)<-1) &....
       (two_D_array2(2,counter_1)>1 | two_D_array2(2,counter_1)<-1)     
    Counter_chnages=Counter_chnages+1;    
    end
    
end

%Write out table of distrubtion of N-termini peptide
fid = fopen('Combinations_of_Tails_and_interactome_which_correlate.txt','w');
fprintf (fid,'%s\t%s\t',...
'Number of combinations which correlate',...
'total Number of combinations'); %Write Header
fprintf(fid,'\n');  
fprintf(fid,'%f	',Counter_chnages(1));
fprintf(fid,'%f	',length(two_D_array2));
fclose(fid);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%check if any of the protein entries matahc the tails data
Comparsion_array=zeros(max(Number_of_terminin_per_proteins),2);   
Observed_termini=zeros(21,2);
Unique_P1_termini=unique(NTAILS.Previousaminoacid);
Sig_Observed_termini_fold_changeA=nan(1,1);
Sig_Observed_termini_fold_changeB=nan(1,1);
Gaus_Observed_termini_fold_change=nan(1,1);
NGaus_Observed_termini_fold_change=nan(1,1);

%create counters
Fold_change_positionA=1;
Fold_change_positionB=1;
NFold_change_positionA=1;
NFold_change_positionB=1;

for max_counter=1:Number_proteins_with_gaussian
    
    %Find how many termini were detected
    [Position_of_accession_num]=ind2sub(size(Proteins_which_overlap(:,1)), strmatch(Proteins_with_Gaussian{max_counter}, Proteins_which_overlap(:,1), 'exact'));
    
    %Find how many termini were detected
    [Gaussian_list_position]=ind2sub(size(Imported_seperated_cyto(:,1)), strmatch(Proteins_with_Gaussian{max_counter}, Imported_seperated_cyto(:,1), 'exact'));
    
    %Check if any of the Gaussian observed for the protein underwhen
    %signficant chnages
    Sign_change_observed=0;
    for subcounter=1:length(Gaussian_list_position)
        
    if isequal(Imported_seperated_cyto(Gaussian_list_position(subcounter),29),{'+'})==1
     Sign_change_observed=1;  
    end
    end
    
    if ~isempty(Position_of_accession_num) & Sign_change_observed==0
    %Find positions of termini equal to value being assessed
    Comparsion_array(Number_of_terminin_per_proteins(Position_of_accession_num),1)=Comparsion_array(Number_of_terminin_per_proteins(Position_of_accession_num),1)+1;
    end
    
    if ~isempty(Position_of_accession_num) & Sign_change_observed==1
    %Find positions of termini equal to value being assessed
    Comparsion_array(Number_of_terminin_per_proteins(Position_of_accession_num),2)=Comparsion_array(Number_of_terminin_per_proteins(Position_of_accession_num),2)+1;
    end
    
    %check what the P1 position is
    if ~isempty(Position_of_accession_num) & Sign_change_observed==0
        
    %Find how many termini were detected
    [Position_in_NTAILS_list]=ind2sub(size(NTAILS.uniprot_accession_num(:,1)), strmatch(Proteins_which_overlap{Position_of_accession_num,2}, NTAILS.uniprot_accession_num(:,1), 'exact'));
    
    for Number_peptide_counter=1:length(Position_in_NTAILS_list)
    for Termini_counter=1:21
        if isequal(Unique_P1_termini(Termini_counter),(NTAILS.Previousaminoacid(Position_in_NTAILS_list(Number_peptide_counter))))
            Observed_termini(Termini_counter,1)=Observed_termini(Termini_counter,1)+1;
        end
    end
    end
    end
    
        %check what the P1 position is
    if ~isempty(Position_of_accession_num) & Sign_change_observed==1
        
    %Find how many termini were detected
    [Position_in_NTAILS_list]=ind2sub(size(NTAILS.uniprot_accession_num(:,1)), strmatch(Proteins_which_overlap{Position_of_accession_num,2}, NTAILS.uniprot_accession_num(:,1), 'exact'));
    
    for Number_peptide_counter=1:length(Position_in_NTAILS_list)
    for Termini_counter=1:21
        if isequal(Unique_P1_termini(Termini_counter),(NTAILS.Previousaminoacid(Position_in_NTAILS_list(Number_peptide_counter))))
            Observed_termini(Termini_counter,2)=Observed_termini(Termini_counter,2)+1;
        end
    end
    end
    end
    
    %check what the average fold change is P1 
    if ~isempty(Position_of_accession_num) & Sign_change_observed==1
        
    %Find how many termini were detected
    [Position_in_NTAILS_list]=ind2sub(size(NTAILS.uniprot_accession_num(:,1)), strmatch(Proteins_which_overlap{Position_of_accession_num,2}, NTAILS.uniprot_accession_num(:,1), 'exact'));
    
    for Number_peptide_counter=1:length(Position_in_NTAILS_list)
    for Termini_counter=1:21
        if isequal(Unique_P1_termini(Termini_counter),(NTAILS.Previousaminoacid(Position_in_NTAILS_list(Number_peptide_counter))))
            Sig_Observed_termini_fold_changeA(Termini_counter,Fold_change_positionA)=str2num(NTAILS.AverageNormalizedMvsL{Position_in_NTAILS_list(Number_peptide_counter)});
            Fold_change_positionA=Fold_change_positionA+1;
            Sig_Observed_termini_fold_changeB(Termini_counter,Fold_change_positionB)=str2num(NTAILS.AverageNormalizedHvsL{Position_in_NTAILS_list(Number_peptide_counter)});
            Fold_change_positionB=Fold_change_positionB+1;
        end
    end
    end
    end
    
    
    %check what the average fold change is P1 
    if ~isempty(Position_of_accession_num) & Sign_change_observed==0
        
    %Find how many termini were detected
    [Position_in_NTAILS_list]=ind2sub(size(NTAILS.uniprot_accession_num(:,1)), strmatch(Proteins_which_overlap{Position_of_accession_num,2}, NTAILS.uniprot_accession_num(:,1), 'exact'));
    
    for Number_peptide_counter=1:length(Position_in_NTAILS_list)
    for Termini_counter=1:21
        if isequal(Unique_P1_termini(Termini_counter),(NTAILS.Previousaminoacid(Position_in_NTAILS_list(Number_peptide_counter))))
            NSig_Observed_termini_fold_changeA(Termini_counter,NFold_change_positionA)=str2num(NTAILS.AverageNormalizedMvsL{Position_in_NTAILS_list(Number_peptide_counter)});
            NFold_change_positionA=NFold_change_positionA+1;
            NSig_Observed_termini_fold_changeB(Termini_counter,NFold_change_positionB)=str2num(NTAILS.AverageNormalizedHvsL{Position_in_NTAILS_list(Number_peptide_counter)});
            NFold_change_positionB=NFold_change_positionB+1;
        end
    end
    end
    end
    
    
      %check what the average fold change is P1 
    if ~isempty(Position_of_accession_num) & Sign_change_observed==1
        
     %Find how many termini were detected
    [Position_of_Gaus]=ind2sub(size(Imported_seperated_cyto(:,1)), strmatch(Proteins_with_Gaussian{max_counter}, Imported_seperated_cyto(:,1), 'exact'));
    
     %Find how many termini were detected
    [Position_in_NTAILS_list]=ind2sub(size(NTAILS.uniprot_accession_num(:,1)), strmatch(Proteins_which_overlap{Position_of_accession_num,2}, NTAILS.uniprot_accession_num(:,1), 'exact'));
    

    for Number_peptide_counterA=1:length(Position_in_NTAILS_list)
     for Number_peptide_counter=1:length(Position_of_Gaus)    
    for Termini_counter=1:21
        if isequal(Unique_P1_termini(Termini_counter),(NTAILS.Previousaminoacid(Position_in_NTAILS_list(Number_peptide_counterA))))
            Gaus_Observed_termini_fold_change(Termini_counter,Fold_change_positionA)=str2num(Imported_seperated_cyto{Position_of_Gaus(Number_peptide_counter),18});
            Fold_change_positionA=Fold_change_positionA+1;
        end
    end
    end
    end
    
    end
    
    
          %check what the average fold change is P1 
    if ~isempty(Position_of_accession_num) & Sign_change_observed==0
        
     %Find how many termini were detected
    [Position_of_Gaus]=ind2sub(size(Imported_seperated_cyto(:,1)), strmatch(Proteins_with_Gaussian{max_counter}, Imported_seperated_cyto(:,1), 'exact'));
    
     %Find how many termini were detected
    [Position_in_NTAILS_list]=ind2sub(size(NTAILS.uniprot_accession_num(:,1)), strmatch(Proteins_which_overlap{Position_of_accession_num,2}, NTAILS.uniprot_accession_num(:,1), 'exact'));
    

    for Number_peptide_counterA=1:length(Position_in_NTAILS_list)
     for Number_peptide_counter=1:length(Position_of_Gaus)    
    for Termini_counter=1:21
        if isequal(Unique_P1_termini(Termini_counter),(NTAILS.Previousaminoacid(Position_in_NTAILS_list(Number_peptide_counterA))))
            NGaus_Observed_termini_fold_change(Termini_counter,Fold_change_positionA)=str2num(Imported_seperated_cyto{Position_of_Gaus(Number_peptide_counter),18});
            Fold_change_positionA=Fold_change_positionA+1;
        end
    end
    end
    end
    
    end
    
    
    
end
  
        
Sig_Observed_termini_fold_changeB(Sig_Observed_termini_fold_changeB==0)=NaN;
Sig_Observed_termini_fold_changeA(Sig_Observed_termini_fold_changeA==0)=NaN;
NSig_Observed_termini_fold_changeB(NSig_Observed_termini_fold_changeB==0)=NaN;
NSig_Observed_termini_fold_changeA(NSig_Observed_termini_fold_changeA==0)=NaN;
Gaus_Observed_termini_fold_change(Gaus_Observed_termini_fold_change==0)=NaN;
NGaus_Observed_termini_fold_change(NGaus_Observed_termini_fold_change==0)=NaN;

%Create X label for figures
Spacing = repmat({' '}, 21, 1);
reshaped_AA=reshape([Unique_P1_termini,Spacing]',[1,42]); 


temp=size(NSig_Observed_termini_fold_changeA);
TempA=reshape(repmat([1:21],length(NSig_Observed_termini_fold_changeA),1)',[1,temp(1)*temp(2)]);
TempB=reshape(NSig_Observed_termini_fold_changeA,[1,temp(1)*temp(2)]);

temp=size(NSig_Observed_termini_fold_changeB);
TempC=reshape(repmat([1:21],length(NSig_Observed_termini_fold_changeB),1)',[1,temp(1)*temp(2)]);
TempD=reshape(NSig_Observed_termini_fold_changeB,[1,temp(1)*temp(2)]);

temp=size(Sig_Observed_termini_fold_changeA);
TempE=reshape(repmat([1:21],length(Sig_Observed_termini_fold_changeA),1)',[1,temp(1)*temp(2)]);
TempF=reshape(Sig_Observed_termini_fold_changeA,[1,temp(1)*temp(2)]);

temp=size(Sig_Observed_termini_fold_changeB);
TempG=reshape(repmat([1:21],length(Sig_Observed_termini_fold_changeB),1)',[1,temp(1)*temp(2)]);
TempH=reshape(Sig_Observed_termini_fold_changeB,[1,temp(1)*temp(2)]);

%Create scatter plot
scatter_plot_1=figure
scatter(TempA,log2(TempB),'filled')
hold on
scatter(TempC,log2(TempD),'filled')
hold on
scatter(TempE,log2(TempF),'filled')
hold on
scatter(TempG,log2(TempH),'filled')
%ylim([-3 3]);
legend('nsig MvsL','nsig HvsL','sig MvsL','sig HvsL')
set(gca,'Xtick',1:21,'XTickLabel',Unique_P1_termini)
print(scatter_plot_1,'-dpdf','-r600', 'Scatterplot_Amino_amino_N-termini_observed_in_proteins_which_chnage_at_gaus_level.pdf');


temp=size(Gaus_Observed_termini_fold_change);
TempA=reshape(repmat([1:21],length(Gaus_Observed_termini_fold_change),1)',[1,temp(1)*temp(2)]);
TempB=reshape(Gaus_Observed_termini_fold_change,[1,temp(1)*temp(2)]);

temp=size(NGaus_Observed_termini_fold_change);
TempC=reshape(repmat([1:21],length(NGaus_Observed_termini_fold_change),1)',[1,temp(1)*temp(2)]);
TempD=reshape(NGaus_Observed_termini_fold_change,[1,temp(1)*temp(2)]);

scatter_plot_2=figure
scatter(TempA,TempB,'filled')
hold on
scatter_plot_1=figure
scatter(TempC,TempD,'filled')
print(scatter_plot_2,'-dpdf','-r600', 'Scatterplot_Amino_amino_N-termini_observed_in_proteins_which_chnage_at_gaus_level2.pdf');

%create normalized array of P1 positions
summed_observed_termini=sum(Observed_termini);
non_sig_Gaus_P1=Observed_termini(:,1)/summed_observed_termini(1);
sig_Gaus_P1=Observed_termini(:,2)/summed_observed_termini(2);

%Create a bar graph
black_array=zeros(1,21);

%Reshape data for figure
B3 = reshape([non_sig_Gaus_P1(:,1),black_array']',[1,42]);
B4 = reshape([black_array',sig_Gaus_P1(:,1)]',[1,42]);

F_signficant=figure
H_B1=bar(B3)
H_B1(1).FaceColor=[colour_to_use(5,:)];
H_B1(1).EdgeColor='k';
hold on
H_B2=bar(B4)
H_B2(1).FaceColor=[colour_to_use(2,:)];
H_B2(1).EdgeColor='k';
title('Amino acid P1 position for N-terminal Acetylated peptides observed to increase');
xlim([0 44])
ylim([0 max(B3)*1.25])

%add X-label
xlabel('Amino acid','Fontsize',10);
set(gca,'Xtick',1:42,'XTickLabel',reshaped_AA)
legend('Proteins shared which do not show changes in Gaussian features','Proteins shared which do show changes in Gaussian features')
print(F_signficant,'-dpdf','-r600', 'Amino_amino_N-termini_observed_in_proteins_which_chnage_at_gaus_level.pdf');


%Write out table of distrubtion of N-termini peptide
fid = fopen('Amino_amino_N-termini_observed_in_proteins_which_chnage_at_gaus_level.txt','w');
fprintf (fid,'%s\t%s\t%s\t',...
'Amino Acid',...
'Proteins shared which do not show changes in Gaussian features',...
'Proteins shared which do show changes in Gaussian features'); %Write Header
fprintf(fid,'\n');
for write_out_counter=1:length(Observed_termini)
fprintf(fid,'%s	',Unique_P1_termini{write_out_counter});    
fprintf(fid,'%f	',Observed_termini(write_out_counter,1));
fprintf(fid,'%f	',Observed_termini(write_out_counter,2));
fprintf(fid,'\n');
end
fclose(fid);





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Create Venn diagram

%Determine the overlap of Interactome to N tails data (All tails data)
I = [length(Unique_proteins_cyto) (Overlap_counter2-1) (Overlap_counter1-1) (length(Unique_N_tails_uniprot)-(Overlap_counter2-1))];
A = [length(Protein_IDs) (length(Unique_proteins_cyto)) (length(Unique_N_tails_uniprot))]; 

F1=figure
h1=venn(A,I,'ErrMinMode','None');
for n=1:length(h1)
set(h1(n),'FaceAlpha', 0.3,'FaceColor',[myC(n,:)],'EdgeColor','black');
end
title('Overlap N-Tails and interactome overlap');
legend('Protein Identified in interactome', 'Protein with protein profiles in interactome','Proteins in N-terminome','FontSize',6, 'Location', 'southoutside')
Output_name5=strcat('Venn_diagram_tails_versus_interactome_overlap.pdf');
print(F1,'-dpdf','-r600', Output_name5);

fid = fopen('Overlap_summary_tails_versus_interactome_overlap.txt','w');
fprintf (fid,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t',...
'Circle A',...
'Circle B',...
'Circle C',...
'Overlap A to B',...
'Overlap A to C',...
'Overlap B to C',...
'Not overlapping in C'); %Write Header
fprintf(fid,'\n');
fprintf(fid,'%f	',A(1));
fprintf(fid,'%f	',A(2));
fprintf(fid,'%f	',A(3));
fprintf(fid,'%f	',I(1));
fprintf(fid,'%f	',I(2));
fprintf(fid,'%f	',I(3));
fprintf(fid,'%f	',I(4));
fclose(fid);


%Show overlap of N-TAILS cytoplasmic interactome and mito interactome
I = [(Overlap_counter2_combined-1) (Overlap_counter2-1) (Overlap_counter2_mito-1) (Overlap_counter3_combined-1)];
A = [length(Protein_IDs) (length(Protein_IDs_mito)) (length(Unique_N_tails_uniprot))]; 

F1=figure
h1=venn(A,I,'ErrMinMode','None');
for n=1:length(h1)
set(h1(n),'FaceAlpha', 0.3,'FaceColor',[myC(n,:)],'EdgeColor','black');
end
title('Overlap Protein identifications N-Tails and interactome overlap');
legend('Protein Identified in Cytosolic interactome', 'Protein Identified in Organelle interactome','Proteins in N-terminome','FontSize',6, 'Location', 'southoutside')
Output_name5=strcat('Venn_diagram_tails_versus_interactomes_cytosol_organelle_overlap.pdf');
print(F1,'-dpdf','-r600', Output_name5);

fid = fopen('Overlap_summary_tails_versus_interactomes_cytosol_organelle.txt','w');
fprintf (fid,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t',...
'Circle A',...
'Circle B',...
'Circle C',...
'Overlap A to B',...
'Overlap A to C',...
'Overlap B to C',...
'Not overlapping in C'); %Write Header
fprintf(fid,'\n');
fprintf(fid,'%f	',A(1));
fprintf(fid,'%f	',A(2));
fprintf(fid,'%f	',A(3));
fprintf(fid,'%f	',I(1));
fprintf(fid,'%f	',I(2));
fprintf(fid,'%f	',I(3));
fprintf(fid,'%f	',I(4));
fclose(fid);


%Determine the overlap of Interactome to N tails data (only quantified)
I2 = [length(Unique_proteins_cyto) length(Overlap_array_for_correlation_NoNAN) (Overlap_counter1_No_NaN-1) (length(Number_unique_Ntails_proteins)-length(Overlap_array_for_correlation_NoNAN-1))];
A2 = [length(Protein_IDs) (length(Unique_proteins_cyto)) (length(Number_unique_Ntails_proteins))]; 

F1=figure
h1=venn(A2,I2,'ErrMinMode','None');
for n=1:length(h1)
set(h1(n),'FaceAlpha', 0.3,'FaceColor',[myC(n,:)],'EdgeColor','black');
end
title('Overlap N-Tails (with Quant measurments) and interactome overlap');
legend('Protein Identified in interactome', 'Protein with protein profiles in interactome','Proteins in N-terminome','FontSize',6, 'Location', 'southoutside')
Output_name5=strcat('Venn_diagram_tails_with_quant_values_versus_interactome_overlap.pdf');
print(F1,'-dpdf','-r600', Output_name5);

fid = fopen('Overlap_summary_tails_versus_interactome_overlap_Only_proteins_with_quant.txt','w');
fprintf (fid,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t',...
'Circle A',...
'Circle B',...
'Circle C',...
'Overlap A to B',...
'Overlap A to C',...
'Overlap B to C',...
'Not overlapping in C'); %Write Header
fprintf(fid,'\n');
fprintf(fid,'%f	',A2(1));
fprintf(fid,'%f	',A2(2));
fprintf(fid,'%f	',A2(3));
fprintf(fid,'%f	',I2(1));
fprintf(fid,'%f	',I2(2));
fprintf(fid,'%f	',I2(3));
fprintf(fid,'%f	',I2(4));
fclose(fid);


%Determine the overlap of Cyto and Mito Interactome to N tails data (only quantified)
I2 = [Overlap_counter1_No_NaN_mito2  (Overlap_counter1_No_NaN-1) (Overlap_counter1_No_NaN_mito-1) Overlap_counter1_Quantfied_all];
A2 = [length(Unique_proteins_cyto) (length(Unique_proteins_mito)) (length(Number_unique_Ntails_proteins))]; 

F1=figure
h1=venn(A2,I2);
for n=1:length(h1)
set(h1(n),'FaceAlpha', 0.3,'FaceColor',[myC(n,:)],'EdgeColor','black');
end
title('Overlap N-Tails (with Quant measurments) and interactomes overlap');
legend('Proteins with protein profiles in Cytosolic interactome','Protein with protein profiles in organelle interactome','Proteins in N-terminome','FontSize',6, 'Location', 'southoutside')
Output_name5=strcat('Venn_diagram_tails_with_quant_values_versus_interactome_mito_cyto_overlap.pdf');
print(F1,'-dpdf','-r600', Output_name5);

fid = fopen('Overlap_summary_tails_versus__mito_cyto_interactome_overlap_Only_proteins_with_quant.txt','w');
fprintf (fid,'%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t',...
'Circle A',...
'Circle B',...
'Circle C',...
'Overlap A to B',...
'Overlap A to C',...
'Overlap B to C',...
'in all samples'); %Write Header
fprintf(fid,'\n');
fprintf(fid,'%f	',A2(1));
fprintf(fid,'%f	',A2(2));
fprintf(fid,'%f	',A2(3));
fprintf(fid,'%f	',I2(1));
fprintf(fid,'%f	',I2(2));
fprintf(fid,'%f	',I2(3));
fprintf(fid,'%f	',I2(4));
fclose(fid);


F1A=figure
venn([length(Unique_proteins_cyto) length(Number_unique_Ntails_proteins) (Overlap_counter1_No_NaN-1)],'FaceColor',{myC(2,:),myC(3,:)},'FaceAlpha',{0.5,0.5},'EdgeColor','black');
title('Overlap N-Tails and all proteins with profiles in interactome overlap');
legend('Protein with protein profiles in interactome', 'Proteins in N-terminome','FontSize',6, 'Location', 'southoutside')
Output_name7=strcat('Venn_diagram_tails_versus_Quantified_Ntails_proteins_vs_Quantified_interactome.pdf');
print(F1A,'-dpdf','-r600', Output_name7);


F2=figure
venn([length(Protein_IDs) (length(Unique_N_tails_uniprot)-(Overlap_counter2-1)) (Overlap_counter2-1)],'FaceColor',{myC(1,:),myC(3,:)},'FaceAlpha',{0.5,0.5},'EdgeColor','black');
title('Overlap N-Tails and all proteins in interactome overlap');
legend('Protein Identified in interactome', 'Proteins in N-terminome','FontSize',6, 'Location', 'southoutside')
Output_name6=strcat('Venn_diagram_tails_versus_all_proteins_in_interactome_overlap.pdf');
print(F2,'-dpdf','-r600', Output_name6);


F3=figure
venn([length(Unique_proteins_cyto) length(Unique_N_tails_uniprot) (Overlap_counter1-1)],'FaceColor',{myC(2,:),myC(3,:)},'FaceAlpha',{0.5,0.5},'EdgeColor','black');
title('Overlap N-Tails and all proteins with profiles in interactome overlap');
legend('Protein with protein profiles in interactome', 'Proteins in N-terminome','FontSize',6, 'Location', 'southoutside')
Output_name7=strcat('Venn_diagram_tails_versus_all_proteins_with_protein_profiles_in_interactome_overlap.pdf');
print(F3,'-dpdf','-r600', Output_name7);



