%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                           Script to Align PCP-SEC-SILAC experiments 
%                                                   Created by Nichollas Scott,
%                                                       Foster lab,UBC, 2014
% Before beginning please place the output from the Gaussian fitting script within the folder containing this script
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Using both Linear and denisty based alignment

%% clear and prep system for use
clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Users defined varibles
User_alignment_window1=20; %for initial alignment with proteins fitted to a single Gaussian
%column based SEC 5 appears to be approperiate, for BN-PAGE determine based
%on results 10
User_alignment_window2=8; %for Second round of alignment
Number_of_experimental_channels=2;  % Defines the number of experiments to be compared
Alignment_to_user= 'MvsL'; %Define the experimental channel to use for alignment       
Image_couter=1; %Start Image counter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Import SILAC ratios
MvsL_filename='Combined_replicates_2014_04_22_contaminates_removed_for_MvsL_scripts.xlsx';
HvsL_filename='Combined_replicates_2014_04_22_contaminates_removed_for_HvsL_scripts.xlsx';
HvsM_filename='Combined_replicates_2014_04_22_contaminates_removed_for_HvsM_scripts.xlsx';
List_of_filename={MvsL_filename,HvsL_filename,HvsM_filename};

%Create Size calibration curve based on Start and end fraction
[SEC_size_alignment]=xlsread('SEC_alignment.xlsx');
SEC_fit=polyfit(SEC_size_alignment(1,:),SEC_size_alignment(2,:),1);

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

%create dictories for script
mkdir('MvsL_alignment');
mkdir('HvsL_alignment');

%import files  
[num_val_MvsL,txt_MvsL] = xlsread(List_of_filename{1}); %Import file raw Maxqaunt output
[num_val_HvsL,txt_HvsL] = xlsread(List_of_filename{2}); %Import file raw Maxqaunt output
[num_val_HvsM,txt_HvsM] = xlsread(List_of_filename{3}); %Import file raw Maxqaunt output

%Number of fractions
fraction_number=size(num_val_MvsL);
fraction_number(2)=fraction_number(2)-1;

% number of fractions,equal to # fractions + 10
Chromatograms_axis=fraction_number(2) +10; 

%lowest value to be shown in Adjusted Chromatograms
User_defined_zero_value=0.2; 
     
Image_couter=1; %Start Image counter

%% Repair any missing values before 
%Replace one missing value with mean of neighbours
%for num_val_MvsL 
for Prot_nr = 1:fraction_number(1) %Defines that proteins should be analyzed, changed to i during testing such be 1:Proteins
    AK1= num_val_MvsL(Prot_nr,:);
    Count =numel(AK1);
    AK2= isnan(AK1); 
    for i=2:(Count-1) 
        if 1==isnan(AK1(i))
                Test1=AK2((i-1):(i+1));
                    if  1==isequal(~Test1,[1 0 1])
                        Number_AK=AK1((i-1):(i+1));
                        Number2_AK1 = Number_AK(~Test1);
                        miss_AK = mean(Number2_AK1); %average if  data point is missing
                        num_val_MvsL(Prot_nr,i) = miss_AK;
                    end
        end
    end

end

%for num_val_HvsL
for Prot_nr = 1:fraction_number(1) %Defines that proteins should be analyzed, changed to i during testing such be 1:Proteins
    AK1= num_val_HvsL(Prot_nr,:);
    AK2= isnan(AK1); 
    for i=2:(Count-1) 
        if 1==isnan(AK1(i))
                Test1=AK2((i-1):(i+1));
                    if  1==isequal(~Test1,[1 0 1])
                        Number_AK=AK1((i-1):(i+1));
                        Number2_AK1 = Number_AK(~Test1);
                        miss_AK = mean(Number2_AK1); %average if  data point is missing
                        num_val_HvsL(Prot_nr,i) = miss_AK;
                    end
        end
    end

end

%for num_val_HvsM
for Prot_nr = 1:fraction_number(1) %Defines that proteins should be analyzed, changed to i during testing such be 1:Proteins
    AK1= num_val_HvsM(Prot_nr,:);
    AK2= isnan(AK1); 
    for i=2:(Count-1) 
        if 1==isnan(AK1(i))
                Test1=AK2((i-1):(i+1));
                    if  1==isequal(~Test1,[1 0 1])
                        Number_AK=AK1((i-1):(i+1));
                        Number2_AK1 = Number_AK(~Test1);
                        miss_AK = mean(Number2_AK1); %average if  data point is missing
                        num_val_HvsM(Prot_nr,i) = miss_AK;
                    end
        end
    end

end


%% Define the experiment to be processed
for Experimental_channel_loop=1:Number_of_experimental_channels
    

    if Experimental_channel_loop ==1
        
        %Create module to divide up files
        %Import Matlab processed data from Gaus script
        Gaus_import_Name= strcat('MvsL\MvsL_Combined_OutputGaus.csv');
        Summary_gausian_infomration_Name=strcat('MvsL\MvsL_Summary_Gausians_for_individual_proteins.csv');
        %import data files Gaus data
        f1 = fopen (Gaus_import_Name);
        Processed_data1 = textscan(f1, '%s', 'Delimiter',',');
        fclose(f1);

        %import data files summary data
        f2 = fopen (Summary_gausian_infomration_Name);
        Processed_data2 = textscan(f2, '%s', 'Delimiter',',');
        fclose(f2);
        
        cd('MvsL_alignment');
        %Create new folders to place output in
        mkdir('Realignment');
        mkdir('Alignment Figures');
        mkdir('Summary tables');
        mkdir('Comparsion tables');
        mkdir('Processed Gaussian');
        mkdir('Combined Realigned data');
        
        Experimental_channel = 'MvsL' % define the MvsL
        
        %Open file to write all image information to
        Image_record_name= strcat('Image Summary.csv');
        Image_record = fopen(Image_record_name,'a+');
        fprintf (Image_record,'%s,%s,\n',...
        'Image number', 'Description');
    
    
        cd('Realignment');
        %Write out table of polynomial fit based on density adjustmetns for replicates
        fid6C = fopen('Summary_table_polynomial_fit_of_density_adjustmetns_for_replicates.csv','a+');
        fprintf (fid6C,'%s,%s,%s,%s,%s\n',...
        'Replicate Number', 'Indeterminates', 'Coefficients', 'R value', 'R squared');  %Write Header
        %Write out table of polynomial fit based on density adjustmetns for replicates
        fid6D = fopen('Summary_table_polynomial_fit_both_adjustmetns_for_replicates.csv','a+');
        fprintf (fid6D,'%s,%s,%s,%s,%s\n',...
        'Replicate Number', 'Indeterminates', 'Coefficients', 'R value', 'R squared');  %Write Header
        cd('../');

    elseif Experimental_channel_loop == 2
        
        %Create module to divide up files
        %Import Matlab processed data from Gaus script
        Gaus_import_Name= strcat('HvsL\HvsL_Combined_OutputGaus.csv');
        Summary_gausian_infomration_Name=strcat('HvsL\HvsL_Summary_Gausians_for_individual_proteins.csv');
        
        %import data files Gaus data
        f1 = fopen (Gaus_import_Name);
        Processed_data1 = textscan(f1, '%s', 'Delimiter',',');
        fclose(f1);

        %import data files summary data
        f2 = fopen (Summary_gausian_infomration_Name);
        Processed_data2 = textscan(f2, '%s', 'Delimiter',',');
        fclose(f2);
        cd('HvsL_alignment');
        mkdir('Realignment');
        mkdir('Alignment Figures');
        mkdir('Summary tables');
        mkdir('Comparsion tables');
        mkdir('Processed Gaussian');
        mkdir('Combined Realigned data');
        
        Experimental_channel = 'HvsL' % define the HvsL
        
        %Open file to write all image information to
        Image_record_name= strcat('Image Summary.csv');
        Image_record = fopen(Image_record_name,'a+');
        fprintf (Image_record,'%s,%s,\n',...
        'Image number', 'Description');
    
        cd('Realignment');
        %Write out table of polynomial fit based on density adjustmetns for replicates
        fid6C = fopen('Summary_table_polynomial_fit_of_density_adjustmetns_for_replicates.csv','a+');
        fprintf (fid6C,'%s,%s,%s,%s,%s\n',...
        'Replicate Number', 'Indeterminates', 'Coefficients', 'R value', 'R squared');  %Write Header

        %Write out table of polynomial fit based on density adjustmetns for replicates
        fid6D = fopen('Summary_table_polynomial_fit_both_adjustmetns_for_replicates.csv','a+');
        fprintf (fid6D,'%s,%s,%s,%s,%s\n',...
        'Replicate Number', 'Indeterminates', 'Coefficients', 'R value', 'R squared');  %Write Header
        cd('../');        
        
    end
        

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

%Reshape data for use in analysis
Summary_gausian_infomration_dimension=size(Processed_data2{:});
Summary_gaussian_infomration=reshape(Processed_data2{:},5,(Summary_gausian_infomration_dimension(1)/5))';

%Determine the number of enteries after reshaped
Gaus_import_correctly_format=size(Gaus_import);
Summary_gausian_correctly_format=size(Summary_gaussian_infomration);

%Determine the number of replicates and unique protein groups
Unique_replicate=length(unique(Gaus_import(2:Gaus_import_correctly_format(1),3)));
Number_unique_proteins=length(unique(Summary_gaussian_infomration(2:Summary_gausian_correctly_format(1),2)));


%Move to folder to output processed Gaussian information into
cd('Processed Gaussian');

%Divid up Combined_OutputGaus_rep
for divider_counter1=1:Unique_replicate
   %Create name of gaussian file to output divided gaus data to
   Process_Gaus_import_Name= strcat(Experimental_channel,'_Combined_OutputGaus_rep',mat2str(divider_counter1),'.csv');
    
   fid_processing= fopen(Process_Gaus_import_Name,'wt'); % create the output file with the header infromation
   fprintf (fid_processing,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n',... %header for OutputGaus output
                'Guassian index number', 'Protein_number','Replicate',...
                'Protein name', 'Height', 'Center', 'Width', 'SSE',...
                'adjrsquare', 'Complex Size'); %Write Header 
    
for divider_counter2=2:Gaus_import_correctly_format(1)
    if str2num(Gaus_import{divider_counter2,3}) ==divider_counter1    
    fprintf(fid_processing,'%6f,%6f,%6f,%s,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,\n',...
         str2num(Gaus_import{divider_counter2,1}),...
         str2num(Gaus_import{divider_counter2,2}),...
         str2num(Gaus_import{divider_counter2,3}),...
         Gaus_import{divider_counter2,4},...
         str2num(Gaus_import{divider_counter2,5}),...
         str2num(Gaus_import{divider_counter2,6}),...
         str2num(Gaus_import{divider_counter2,7}),...
         str2num(Gaus_import{divider_counter2,8}),...
         str2num(Gaus_import{divider_counter2,9}),...
         str2num(Gaus_import{divider_counter2,10}));
    end
end
fclose(fid_processing);
end

          
%Divid up Summary_Gausians_for_individual_proteins_rep
for divider_counter1=1:Unique_replicate
   %Create name of gaussian file to output divided gaus data to
   Process_Summary_gausian_info_Name= strcat(Experimental_channel,'_Summary_Gausians_for_individual_proteins_rep',mat2str(divider_counter1),'.csv');
   
   fid_processing2= fopen(Process_Summary_gausian_info_Name,'wt'); % create the output file with the header infromation
   fprintf (fid_processing2,'%s,%s,%s,%s,%s\n',... %header for OutputGaus output
                'Protein_number', 'Protein_name','Number_of_Gausians_detected',...
                'Number_of_Gausians_within_defined_boundaries',...
                'Number_of_Gausians_filtered'); %Write Header 
    				
for divider_counter2=2+(Number_unique_proteins*(divider_counter1-1)):(Number_unique_proteins+(Number_unique_proteins*(divider_counter1-1))+1)
      fprintf(fid_processing2,'%6f,%s,%6f,%6f,%6f,\n',...
      str2num(Summary_gaussian_infomration{divider_counter2,1}),...
      Summary_gaussian_infomration{divider_counter2,2},...
      str2num(Summary_gaussian_infomration{divider_counter2,3}),...
      str2num(Summary_gaussian_infomration{divider_counter2,4}),...
      str2num(Summary_gaussian_infomration{divider_counter2,5}));
end

fclose(fid_processing2);   
end
       
%Divid up Summary_Gausians_for_individual_proteins_rep
for divider_counter1=1:Unique_replicate
    
    %Create name of gaussian file to output divided gaus data for MvsL
    Process_Replicate_raw_data1= strcat('MvsL_Raw_data_maxquant_rep',mat2str(divider_counter1),'.csv');
   
    fid_processing3= fopen(Process_Replicate_raw_data1,'wt'); % create the output file with the header infromation
    fprintf(fid_processing3,'%s, ', txt_MvsL{1,1:fraction_number(2)+2}); %header for OutputGaus output
   	fprintf(fid_processing3,'\n');			
    for divider_counter3=1+(Number_unique_proteins*(divider_counter1-1)):(Number_unique_proteins+(Number_unique_proteins*(divider_counter1-1)))
    fprintf(fid_processing3,'%s, %6.4f,', txt_MvsL{divider_counter3+1,1}, num_val_MvsL(divider_counter3,1));
    fprintf(fid_processing3,'%6.4g,',num_val_MvsL(divider_counter3,(2:fraction_number(2)+1))); %Chromatogram information
    fprintf(fid_processing3,'\n');
    end
fclose(fid_processing3); 

    %Create name of gaussian file to output divided gaus data for HvsL
    Process_Replicate_raw_data2= strcat('HvsL_Raw_data_maxquant_rep',mat2str(divider_counter1),'.csv');
   
    fid_processing4= fopen(Process_Replicate_raw_data2,'wt'); % create the output file with the header infromation
    fprintf (fid_processing4,'%s, ', txt_HvsL{1,1:fraction_number(2)+2}); %header for OutputGaus output
   	fprintf(fid_processing4,'\n');			
    for divider_counter3=1+(Number_unique_proteins*(divider_counter1-1)):(Number_unique_proteins+(Number_unique_proteins*(divider_counter1-1)))
    fprintf(fid_processing4,'%s, %6.4f,', txt_HvsL{divider_counter3+1,1}, num_val_HvsL(divider_counter3,1));
    fprintf(fid_processing4,'%6.4g,',num_val_HvsL(divider_counter3,(2:fraction_number(2)+1))); %Chromatogram information
    fprintf(fid_processing4,'\n');
    end
fclose(fid_processing4); 

    %Create name of gaussian file to output divided gaus data for HvsL
    Process_Replicate_raw_data3= strcat('HvsM_Raw_data_maxquant_rep',mat2str(divider_counter1),'.csv');
   
    fid_processing5= fopen(Process_Replicate_raw_data3,'wt'); % create the output file with the header infromation
    fprintf (fid_processing5,'%s, ', txt_HvsM{1,1:fraction_number(2)+2}); %header for OutputGaus output
   	fprintf(fid_processing5,'\n');			
    for divider_counter3=1+(Number_unique_proteins*(divider_counter1-1)):(Number_unique_proteins+(Number_unique_proteins*(divider_counter1-1)))
    fprintf(fid_processing5,'%s, %6.4f,', txt_HvsM{divider_counter3+1,1}, num_val_HvsM(divider_counter3,1));
    fprintf(fid_processing5,'%6.4g,',num_val_HvsM(divider_counter3,(2:fraction_number(2)+1))); %Chromatogram information
    fprintf(fid_processing5,'\n');
    end
fclose(fid_processing5); 

end

cd('../');

%%%%%%%%%%%%%%%%%%%%%%%
%% Start loop to align and correct data
for replicates= 1:Unique_replicate
    
%Import Matlab processed data
Gaus_import_Name= strcat('Processed Gaussian\',Experimental_channel,'_Combined_OutputGaus_rep',mat2str(replicates),'.csv');
Gaus_import=importdata(Gaus_import_Name);           %Import the Gaus data from previous scripts
Summary_gausian_infomration_Name=strcat('Processed Gaussian\',Experimental_channel,'_Summary_Gausians_for_individual_proteins_rep',mat2str(replicates),'.csv');
Summary_gausian_infomration=importdata(Summary_gausian_infomration_Name); %Import summary data from previous scripts

%Write out varibles
summerised_protein_number=cell(1,1);
summerised_names=cell(1,1);
summerised_guas_number=zeros(1,1);

%Identify proteins with gausian detected, Although the filter script does this
%this will remove proteins previously determined to not have gaussian
%within the boundaries set by the user.
check=length(Summary_gausian_infomration.data);
checker_counter=1;
for checker= 1:check
    values=Summary_gausian_infomration.data(checker,2);
    values=values>0;
    if values==1;
        summerised_protein_number(checker_counter,1)=Summary_gausian_infomration.textdata(checker+1,1);
        summerised_names(checker_counter,1)=Summary_gausian_infomration.textdata(checker+1,2);
        summerised_guas_number(checker_counter,1)=Summary_gausian_infomration.data(checker,2);
        checker_counter=checker_counter+1;
    end
end

%Write out file of proteins with gaussians within boundaries set by the user.
fid_Name= strcat('Summary_Gausians_for_individual_proteins_with_guassiane_rep',mat2str(replicates),'.csv');
fid = fopen(fid_Name,'w');
fprintf (fid,'%s,%s,%s\n',...
'Protein_number', 'Gene_name', 'Number_of_Gausians_detected');               %Write Header
for i=1:(checker_counter-1)
    fprintf(fid, '%s,%s,%6.4f\n',...
        summerised_protein_number{i},summerised_names{i},summerised_guas_number(i));
end
fclose(fid);

%Move output to folder
movefile(fid_Name, 'Summary tables');

%Set varibles to zero
summerised_protein_number_G1=cell(1,1);
summerised_names_G1=cell(1,1);
summerised_guas_number_G1=zeros(1,1);

%Identify proteins with only one gausian detected
checker_counter2=1;
for checker= 1:check
    values=Summary_gausian_infomration.data(checker,2);
    if 1== values; %%%there seems to be a problem here
        summerised_protein_number_G1(checker_counter2,1)=Summary_gausian_infomration.textdata(checker+1,1);
        summerised_names_G1(checker_counter2,1)=Summary_gausian_infomration.textdata(checker+1,2);
        summerised_guas_number_G1(checker_counter2,1)=Summary_gausian_infomration.data(checker,2);
        checker_counter2=checker_counter2+1;
    end
end

%Write out file of proteins with gaussians within boundaries set by the user.
fid2_Name= strcat('Summary_Gausians_for_individual_proteins_with_one_guassiane_rep',mat2str(replicates),'.csv');
fid2 = fopen(fid2_Name,'w');
fprintf (fid2,'%s,%s,%s\n',...
'Protein_number', 'Gene_name', 'Number_of_Gausians_detected');               %Write Header
for i=1:(checker_counter2-1)
    fprintf(fid2, '%s,%s,%6.4f\n',...
        summerised_protein_number_G1{i},summerised_names_G1{i},summerised_guas_number_G1(i));
end
fclose(fid2);

%Move output to folder
movefile(fid2_Name, 'Summary tables');

%Create varibles to write to
Table_value_for_proteins_one_gaus=zeros(1,5);

%Make table of all the values from proteins with only single gaussian
check2=length(summerised_names_G1);
Table_value_for_proteins_one_gaus=zeros(check2,5);
for ii=1:check2
    find_string= strmatch(summerised_names_G1{ii},Gaus_import.textdata(:,4), 'exact'); 
    Table_value_for_proteins_one_gaus(ii,1)=Gaus_import.data(find_string-1,1); %minus 1 as .data does not contain the header row found in .textdata
    Table_value_for_proteins_one_gaus(ii,2)=Gaus_import.data(find_string-1,2); %minus 1 as .data does not contain the header row found in .textdata
    Table_value_for_proteins_one_gaus(ii,3)=Gaus_import.data(find_string-1,3); %minus 1 as .data does not contain the header row found in .textdata
    Table_value_for_proteins_one_gaus(ii,4)=Gaus_import.data(find_string-1,4); %minus 1 as .data does not contain the header row found in .textdata
    Table_value_for_proteins_one_gaus(ii,5)=Gaus_import.data(find_string-1,5); %minus 1 as .data does not contain the header row found in .textdata
end
   
%Clear structure before write
clear summerised_names_G1_values;
clear Sorted_structure;

%Set up a structure array of protein with single guassian
for iv=1:check2
    summerised_names_G1_values(iv).Protein_name = summerised_names_G1{iv};
    summerised_names_G1_values(iv).Height = Table_value_for_proteins_one_gaus(iv,1);
    summerised_names_G1_values(iv).Center = Table_value_for_proteins_one_gaus(iv,2);
    summerised_names_G1_values(iv).Width = Table_value_for_proteins_one_gaus(iv,3);
    summerised_names_G1_values(iv).SSE = Table_value_for_proteins_one_gaus(iv,4);
    summerised_names_G1_values(iv).adjrsquare = Table_value_for_proteins_one_gaus(iv,5);
end

%sort the structure array
[location,location]=sort([summerised_names_G1_values.Center]);
Sorted_structure=summerised_names_G1_values(location);

%write out sorted structure array
fid3_Name = strcat('Sorted_Gausians_for_individual_proteins_with_one_guassiane_rep',mat2str(replicates),'_','.csv');
fid3 = fopen(fid3_Name,'w');
fprintf (fid3,'%s,%s,%s,%s,%s,%s\n',...
'Protein_name', 'Height', 'Center','Width','SSE','adjrsquare');               %Write Header
for v=1:check2
    fprintf(fid3, '%s,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,\n',...
        Sorted_structure(v).Protein_name, Sorted_structure(v).Height,...
        Sorted_structure(v).Center, Sorted_structure(v).Width,...
        Sorted_structure(v).SSE, Sorted_structure(v).adjrsquare);
end
fclose(fid3);

%Move output to folder
movefile(fid3_Name, 'Summary tables');

end %close loop for replicates

%%Identify protiens with only a signle gaussain detected in each experiment
Summary_details_replicate.intersect=zeros(Unique_replicate,Unique_replicate);

for replicates= 1:Unique_replicate
    
%impott list of proteins with a single guassian from replicates to test
cd('Summary tables');
Single_Gaus_import_Name= strcat('Sorted_Gausians_for_individual_proteins_with_one_guassiane_rep',mat2str(replicates),'_.csv');
Single_guassian_import=importdata(Single_Gaus_import_Name);
number_of_protein_one_gaus=length(Single_guassian_import.textdata);
Protein_names_gaussian=Single_guassian_import.textdata(2:number_of_protein_one_gaus);
cd('../');

for replicates2= 1:Unique_replicate
%impott list of proteins to compared to list abover
cd('Summary tables');
Single_Gaus_import_Name_replicate_to_test= strcat('Sorted_Gausians_for_individual_proteins_with_one_guassiane_rep',mat2str(replicates2),'_.csv');
Single_guassian_import_to_compare=importdata(Single_Gaus_import_Name_replicate_to_test);
number_of_protein_one_gaus_to_compare=length(Single_guassian_import_to_compare.textdata);
Protein_names_gaussian_to_compare=Single_guassian_import_to_compare.textdata(2:number_of_protein_one_gaus_to_compare);
cd('../');

%Compare replicates to determine the best replicate to align other
%replicates to
Summary_details_replicate.intersect(replicates(1),replicates2(1))=length(intersect(Protein_names_gaussian,Protein_names_gaussian_to_compare));
Summary_details_replicate.intersect(logical(eye(size(Summary_details_replicate.intersect)))) = 0 % remove the diagonal values
end
end

%Determine which replicate in intersect summary is best
for i=1:Unique_replicate
    maximum_intersects(i)=sum(Summary_details_replicate.intersect(i,:));
end

%Find replicate with the maximum number of gaussians found in other
%replicates
[num idx] = max(maximum_intersects(:));
[x_maximum y_maximum] = ind2sub(size(maximum_intersects),idx);

%remove the replicate that will be used for the alignment
replicates_minus_alignment=[1:1:Unique_replicate]
replicates_minus_alignment(y_maximum)=[];

%Write out which replicate were alignment against each other
cd('Realignment');
fid_Alignment_Name= fopen('Alignment_table.csv','w');
fprintf (fid_Alignment_Name,'%s,%s\n',...
'Replicate adjusted', 'Aligned against replicate');  %Write Header
for fid_Alignment_Name_counter=1:Unique_replicate-1
    fprintf(fid_Alignment_Name, '%6.4f,%6.4f\n',...
        replicates_minus_alignment(fid_Alignment_Name_counter), y_maximum);
end
fclose(fid_Alignment_Name);
cd('../');

%Create matrix to deposit the calculated correction values (based on
%polyfit)
Polyfit_summary_user=zeros(Unique_replicate-1,2);

%y_maximum is the replicate which the other replicates are aligned to
for replicate_to_align_against= y_maximum
    
%impott list of proteins with a single guassian from replicates to test
cd('Summary tables');
Single_Gaus_import_Name= strcat('Sorted_Gausians_for_individual_proteins_with_one_guassiane_rep',mat2str(replicate_to_align_against),'_.csv');
Single_guassian_import=importdata(Single_Gaus_import_Name);
number_of_protein_one_gaus=length(Single_guassian_import.textdata);
Protein_names_gaussian=Single_guassian_import.textdata(2:number_of_protein_one_gaus);
cd('../');

for replicates1= 1:Unique_replicate-1
%Call the replicates to check
replicates2=replicates_minus_alignment(replicates1);

%impott list of proteins to compared to list abover
cd('Summary tables');
Single_Gaus_import_Name_replicate_to_test= strcat('Sorted_Gausians_for_individual_proteins_with_one_guassiane_rep',mat2str(replicates2),'_.csv');
Single_guassian_import_to_compare=importdata(Single_Gaus_import_Name_replicate_to_test);
number_of_protein_one_gaus_to_compare=length(Single_guassian_import_to_compare.textdata);
Protein_names_gaussian_to_compare=Single_guassian_import_to_compare.textdata(2:number_of_protein_one_gaus_to_compare);
cd('../');


%combine list and calculate difference in values
Combined_list=intersect(Protein_names_gaussian,Protein_names_gaussian_to_compare);
Combined_list_lenght=length(Combined_list);
for iv =1:Combined_list_lenght
    protein_name_in_both_lists=Combined_list{iv};
    [location_rep1]=ind2sub(size(Protein_names_gaussian), strmatch(protein_name_in_both_lists, Protein_names_gaussian, 'exact'));
    [location_rep2]=ind2sub(size(Protein_names_gaussian_to_compare), strmatch(protein_name_in_both_lists, Protein_names_gaussian_to_compare, 'exact'));
    %copy values from replicate one list
     Table_value_for_rep1(iv,1)=Single_guassian_import.data(location_rep1,1);
     Table_value_for_rep1(iv,2)=Single_guassian_import.data(location_rep1,2);
     Table_value_for_rep1(iv,3)=Single_guassian_import.data(location_rep1,3);
     Table_value_for_rep1(iv,4)=Single_guassian_import.data(location_rep1,4);
     Table_value_for_rep1(iv,5)=Single_guassian_import.data(location_rep1,5);
    %copy values from replicate two list
     Table_value_for_rep2(iv,1)=Single_guassian_import_to_compare.data(location_rep2,1);
     Table_value_for_rep2(iv,2)=Single_guassian_import_to_compare.data(location_rep2,2);
     Table_value_for_rep2(iv,3)=Single_guassian_import_to_compare.data(location_rep2,3);
     Table_value_for_rep2(iv,4)=Single_guassian_import_to_compare.data(location_rep2,4);
     Table_value_for_rep2(iv,5)=Single_guassian_import_to_compare.data(location_rep2,5);
    %Difference between height, center and width values
     Delta_height(iv,1)= Table_value_for_rep2(iv,1)-Table_value_for_rep1(iv,1);
     Delta_center(iv,1)= Table_value_for_rep2(iv,2)-Table_value_for_rep1(iv,2);
     Delta_width(iv,1) = Table_value_for_rep2(iv,3)-Table_value_for_rep1(iv,3);
    for Chromatograms_axis_counter =1:Chromatograms_axis
    Adjusted_Chromatograms_rep1(Chromatograms_axis_counter)... %Create gaussian curve from fitted data
        =  Table_value_for_rep1(iv,1)*exp(-((Chromatograms_axis_counter-(Table_value_for_rep1(iv,2)+5))/Table_value_for_rep1(iv,3)).^2);
    Adjusted_Chromatograms_rep2(Chromatograms_axis_counter)... %Create gaussian curve from fitted data
        =  Table_value_for_rep2(iv,1)*exp(-((Chromatograms_axis_counter-(Table_value_for_rep2(iv,2)+5))/Table_value_for_rep2(iv,3)).^2);
    end
    Adjusted_Chromatograms_rep1(Adjusted_Chromatograms_rep1<0.1) = User_defined_zero_value(1);
    Adjusted_Chromatograms_rep2(Adjusted_Chromatograms_rep2<0.1) = User_defined_zero_value(1);
    EuDis(iv,1) = sqrt(sum((Adjusted_Chromatograms_rep1 - Adjusted_Chromatograms_rep2) .^ 2));
end

%Write out data files
fid4_Name= strcat('Comparsion of guassian properties between replicates_rep',mat2str(replicate_to_align_against),'_compared to_rep',mat2str(replicates2),'.csv')
fid4 = fopen(fid4_Name,'w');
fprintf (fid4,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n',...
'Protein_name', 'Height repliacte 1', 'Center repliacte 1','Width repliacte 1',...
'SSE repliacte 1','adjrsquare repliacte 1','Height repliacte 2',... 
'Center repliacte 2','Width repliacte 2', 'SSE repliacte 2',...
'adjrsquare repliacte 2','Delta height', 'Delta center','Delta width','Euclidean distance');  %Write Header
for iiv=1:Combined_list_lenght
    fprintf(fid4, '%s,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f\n',...
        Combined_list{iiv}, Table_value_for_rep1(iiv,1),...
        Table_value_for_rep1(iiv,2), Table_value_for_rep1(iiv,3),...
        Table_value_for_rep1(iiv,4), Table_value_for_rep1(iiv,5),...
        Table_value_for_rep2(iiv,1), Table_value_for_rep2(iiv,2),...
        Table_value_for_rep2(iiv,3), Table_value_for_rep2(iiv,4),...
        Table_value_for_rep2(iiv,5), Delta_height(iiv,1),...
        Delta_center(iiv,1), Delta_width(iiv,1), EuDis(iiv,1));
end
fclose(fid4);

%Move output to comparsion folder
movefile(fid4_Name, 'Comparsion tables');

%Sort proteins with one guassian in both replicates to enable the
%comparsion of the center
linear_counter1=1;
for iiiv =1:Combined_list_lenght
    sort_array(linear_counter1).Protein_name = Combined_list{iiiv};
    sort_array(linear_counter1).Height_rep1 = Table_value_for_rep1(iiiv,1);
    sort_array(linear_counter1).Height_rep2 = Table_value_for_rep2(iiiv,1);
    sort_array(linear_counter1).Center_rep1 = Table_value_for_rep1(iiiv,2);
    sort_array(linear_counter1).Center_rep2 = Table_value_for_rep2(iiiv,2);
    sort_array(linear_counter1).Width_rep1 = Table_value_for_rep1(iiiv,3);
    sort_array(linear_counter1).Width_rep2 = Table_value_for_rep2(iiiv,3);
    sort_array(linear_counter1).Delta_center = Delta_center(iiiv,1);
    sort_array(linear_counter1).Delta_height = Delta_height(iiiv,1);
    sort_array(linear_counter1).Delta_width = Delta_width(iiiv,1);
    sort_array(linear_counter1).Euclidean_distance = EuDis(iiiv,1);
    linear_counter1=linear_counter1+1;
end

%Sort structured array by replicate 1 center
[location2,location2]=sort([sort_array.Center_rep1]);
Sorted_protein_guassian_array=sort_array(location2);
Sorted_protein_guassian_array2=sort_array(location2);

%reset rows to remove
Rows_to_remove=zeros(1,1);

for Counter_ten = 1:Combined_list_lenght
if Sorted_protein_guassian_array(Counter_ten).Center_rep1 > fraction_number(2) | Sorted_protein_guassian_array(Counter_ten).Center_rep2 > fraction_number(2)
Rows_to_remove(Counter_ten)=1;
end
end

%Remove rows wiuth values over user define threshold
Sorted_protein_guassian_array(Rows_to_remove(:)==1)=[];

%Determine lenght of array with centers less the user defined threshold
Combined_list_lenght_user=length(Sorted_protein_guassian_array);
Combined_list_lenght_user1=length(Sorted_protein_guassian_array2);

%%Manipulate data to use for figure of the alginment of protein with single gaussians
for Counter_ten = 1:Combined_list_lenght_user
Protein_name_within_tolerance1{Counter_ten}= Sorted_protein_guassian_array(Counter_ten).Protein_name; %copy data from structure array to double array    
x_values1(Counter_ten)=Sorted_protein_guassian_array(Counter_ten).Center_rep1; %copy data from structure array to double array
y_values1(Counter_ten)=Sorted_protein_guassian_array(Counter_ten).Center_rep2; %copy data from structure array to double array
z_values1(Counter_ten)=Sorted_protein_guassian_array(Counter_ten).Delta_center;%copy data from structure array to double array
end

for Counter_ten1 = 1:Combined_list_lenght_user1
x_values2(Counter_ten1)=Sorted_protein_guassian_array2(Counter_ten1).Center_rep1; %copy data from structure array to double array
y_values2(Counter_ten1)=Sorted_protein_guassian_array2(Counter_ten1).Center_rep2; %copy data from structure array to double array
z_values2(Counter_ten1)=Sorted_protein_guassian_array2(Counter_ten1).Delta_center;%copy data from structure array to double array
end

%calculate Polynomial curve that fits data
Curve_user_defined1=polyfit(x_values1,y_values1,1);
polyfit_x_value1=round(Curve_user_defined1(1)*100000)/100000;
polyfit_consent1=round(Curve_user_defined1(2)*100000)/100000;
Equation_of_polyfit1=strcat('Equation Of The Line: ', mat2str(polyfit_x_value1),'x +',mat2str(polyfit_consent1));
%testing Correlation Coefficients
correlation_window=corrcoef(x_values1,y_values1);
R_value1=round(correlation_window(2)*100000)/100000;
R_squared1=round((correlation_window(2)*correlation_window(2))*1000)/1000;
R_value_for_figure1=strcat('Correlation Coefficient, R : ', mat2str(R_value1));
R_squared_for_figure1=strcat('Coefficient of Determination, R^2 : ', mat2str(R_squared1));

%Save polyfit results to file to write out
Polyfit_summary_global(replicates1,1)=polyfit_x_value1;
Polyfit_summary_global(replicates1,2)=polyfit_consent1;
R_values_summary_global(replicates1,1)=R_value1;
R_values_summary_global(replicates1,2)=R_squared1;

%%images for all proteins shared between replicates
%3D scatter plot showing the delta difference between replicates
Combined_delta=log10((abs(Delta_height))+(abs(Delta_width))+abs(Delta_center)); %measure of total divergence between replicates
%Note the delta heigth has been weigthed (the *5) to account for the
%magitude of difference with this value compared to delta width and delta center 
image_3D_scatter_pr=figure;
image=scatter3(abs(Delta_center), abs(Delta_height), abs(Delta_width), 50, Combined_delta,'filled');
colourbar_figure1= colorbar('location', 'SouthOutside');
xlabel(colourbar_figure1, 'Divergence');
% Add title and axis labels
title('Properties of Gaussian between replicates');
xlabel('Delta Center');
ylabel('Delta Height');
zlabel('Delta Width');
view(-45, 20);
ThreeD_scatter_plot_name=strcat('Comparsion of the one guassian properties between replicate',mat2str(replicate_to_align_against),' compared to replicate',mat2str(replicates2));

%generate image name
Image_name=strcat('image_',mat2str(Image_couter),'.png');

%Save image
saveas(image, Image_name);

%Save file name to image record
fprintf (Image_record,'%s,%s,\n',...
Image_name, ThreeD_scatter_plot_name);

%Add one to image counter
Image_couter=Image_couter+1;

%Move output to Figure folder
movefile(Image_name, 'Alignment Figures');

%Figure of correlation between centers of raw data
Alignment_graphic= figure;
scatter(x_values2, y_values2, 50,  abs(z_values2), 'filled');
lsline
title({'Comparsion of proteins with single guassian between replicates'});
xlabel_name= strcat('Biological replicate',mat2str(replicate_to_align_against));
ylabel_name= strcat('Biological replicate',mat2str(replicates2));
xlabel(xlabel_name);
ylabel(ylabel_name);
%Add text with R, R2 and equation of the line
text(7, 48, Equation_of_polyfit1);
text(7, 45, R_value_for_figure1);
text(7, 42, R_squared_for_figure1);
%Add a colorbar with tick labels
colourbar_figure= colorbar('location', 'EastOutside');
ylabel(colourbar_figure, 'Delta Center');
line_plot_name=strcat('Comparsion of all proteins with single guassian between replicate',mat2str(replicate_to_align_against),'to',mat2str(replicates2));

%generate image name
Image_name=strcat('image_',mat2str(Image_couter),'.png');

%Save image
saveas(Alignment_graphic, Image_name);

%Save file name to image record
fprintf (Image_record,'%s,%s,\n',...
Image_name, line_plot_name);

%Add one to image counter
Image_couter=Image_couter+1;

%Move output to Figure folder
movefile(Image_name, 'Alignment Figures');


%Histogram of divergence, enables rapid comparsion of how similar
%replicates are
fh = figure;
raw_divergent=(abs(Delta_height))+(abs(Delta_width))+abs(Delta_center);
raw_bin_values=1:50;
hist(raw_divergent,raw_bin_values)
title('Divergence between replicates');
xlabel('Divergence');
ylabel('Number of identifications');
histo_plot_name=strcat('Histogram of divergence between replicate',mat2str(replicate_to_align_against),'_compared to_replicate',mat2str(replicates2));

%generate image name
Image_name=strcat('image_',mat2str(Image_couter),'.png');

%Save image
saveas(fh, Image_name, 'png');
close(fh)

%Save file name to image record
fprintf (Image_record,'%s,%s,\n',...
Image_name, histo_plot_name);

%Add one to image counter
Image_couter=Image_couter+1;

%Move output to Figure folder
movefile(Image_name, 'Alignment Figures');


%%Check linearity, check the equ of the line. 
%the height and center should be linked
linearity_check = User_alignment_window1; % user defined setting of the time window to use
Linearity_check_Delta_center=Delta_center;
Linearity_check_Delta_center(Linearity_check_Delta_center>linearity_check)=NaN;
Linearity_check_Delta_center(Linearity_check_Delta_center<-linearity_check)=NaN;
Logic_Linearity_check_Delta_center=~isnan(Linearity_check_Delta_center(:,1));
linear_counter2=1;
for iiiv =1:Combined_list_lenght
    if Logic_Linearity_check_Delta_center(iiiv) == 1
    linear_check_array(linear_counter2).Protein_name = Combined_list{iiiv};
    linear_check_array(linear_counter2).Height_rep1 = Table_value_for_rep1(iiiv,1);
    linear_check_array(linear_counter2).Height_rep2 = Table_value_for_rep2(iiiv,1);
    linear_check_array(linear_counter2).Center_rep1 = Table_value_for_rep1(iiiv,2);
    linear_check_array(linear_counter2).Center_rep2 = Table_value_for_rep2(iiiv,2);
    linear_check_array(linear_counter2).Width_rep1 = Table_value_for_rep1(iiiv,3);
    linear_check_array(linear_counter2).Width_rep2 = Table_value_for_rep2(iiiv,3);
    linear_check_array(linear_counter2).Delta_center = Delta_center(iiiv,1);
    linear_check_array(linear_counter2).Delta_height = Delta_height(iiiv,1);
    linear_check_array(linear_counter2).Delta_width = Delta_width(iiiv,1);
    linear_check_array(linear_counter2).Euclidean_distance = EuDis(iiiv,1);    
    linear_counter2=linear_counter2+1;
    end
end

%Sort structured array by replicate 1 center
[location3,location3]=sort([linear_check_array.Center_rep1]);
Sorted_linearity_check=linear_check_array(location3);
Sorted_linearity_check2=linear_check_array(location3);

%reset rows to remove
Rows_to_remove=zeros(1,1);

for Counter_ten = 1:linear_counter2-1
if Sorted_linearity_check(Counter_ten).Center_rep1 > fraction_number(2) | Sorted_linearity_check(Counter_ten).Center_rep2 > fraction_number(2)
Rows_to_remove(Counter_ten)=1;
end
end

%Remove rows wiuth values over user define threshold
Sorted_linearity_check(Rows_to_remove(:)==1)=[];

%Determine lenght of array with centers less the user defined threshold
Linear_list_lenght_user=length(Sorted_linearity_check);
Linear_list_lenght_user1=length(Sorted_linearity_check2);

%%Create figure of the alginment of protein with single gaussians and
%%within user defined tolerance window (inital set to 5)
%copy data to use in plot
for Counter_ten = 1:Linear_list_lenght_user
Protein_name_within_tolerance{Counter_ten}= Sorted_linearity_check(Counter_ten).Protein_name; %copy data from structure array to double array    
x_values3(Counter_ten)=Sorted_linearity_check(Counter_ten).Center_rep1; %copy data from structure array to double array
y_values3(Counter_ten)=Sorted_linearity_check(Counter_ten).Center_rep2; %copy data from structure array to double array
z_values3(Counter_ten)=Sorted_linearity_check(Counter_ten).Delta_center;%copy data from structure array to double array
Delta_height2(Counter_ten)=Sorted_linearity_check(Counter_ten).Delta_height;
Delta_width2(Counter_ten)=Sorted_linearity_check(Counter_ten).Delta_width;
Delta_center2(Counter_ten)=Sorted_linearity_check(Counter_ten).Delta_center;
Euclidean_distance2(Counter_ten)=Sorted_linearity_check(Counter_ten).Euclidean_distance;
end

%%Create arrays to use for figure figure
for Counter_ten = 1:Linear_list_lenght_user1   
x_values4(Counter_ten)=Sorted_linearity_check2(Counter_ten).Center_rep1; %copy data from structure array to double array
y_values4(Counter_ten)=Sorted_linearity_check2(Counter_ten).Center_rep2; %copy data from structure array to double array
z_values4(Counter_ten)=Sorted_linearity_check2(Counter_ten).Delta_center;%copy data from structure array to double array
end

%Sorting of Delta values to determine line for precentile
%height
Delta_height3=sort(abs(Delta_height2));

%width
Delta_width3=sort(abs(Delta_width2));

%center
Delta_center3=sort(abs(Delta_center2));

%Euclidean distance
Euclidean_distance3=sort(Euclidean_distance2);


%Write out table of proteins with a single Gaussian in both replicates with
%Center within user defined tolerance
cd('Comparsion tables');
fid5_name=strcat('Comparsion_of_Single_guassian_between_replicates',mat2str(replicate_to_align_against),'to',mat2str(replicates2),'_data_within_user_settings.csv');
fid5 = fopen(fid5_name,'w');
fprintf (fid5,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n',...
'Protein_name', 'Center repliacte 1', 'Center repliacte 2','Delta center',...
'Height repliacte 1','Height repliacte 2', 'Delta height','Width repliacte 1',...
'Width repliacte 2', 'Delta Width', 'Euclidean distance');  %Write Header
for ix=1:Linear_list_lenght_user
    fprintf(fid5, '%s,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f,%6.4f\n',...
        Protein_name_within_tolerance{ix}, Sorted_linearity_check(ix).Center_rep1,...
        Sorted_linearity_check(ix).Center_rep2, Sorted_linearity_check(ix).Delta_center,...
        Sorted_linearity_check(ix).Height_rep1, Sorted_linearity_check(ix).Height_rep2,...
        Sorted_linearity_check(ix).Delta_height, Sorted_linearity_check(ix).Width_rep1,...
        Sorted_linearity_check(ix).Width_rep2, Sorted_linearity_check(ix).Delta_width,...
        Sorted_linearity_check(ix).Euclidean_distance);
end
fclose(fid5);
cd('../');

%calculate Polynomial curve that fits data
Curve_user_defined=polyfit(x_values3,y_values3,1);
polyfit_x_value=round(Curve_user_defined(1)*100000)/100000;
polyfit_consent=round(Curve_user_defined(2)*100000)/100000;
Equation_of_polyfit=strcat('Equation Of The Line: ', mat2str(polyfit_x_value),'x ',mat2str(polyfit_consent));
%testing Correlation Coefficients
correlation_user_defined_window=corrcoef(x_values3,y_values3);
R_value=round(correlation_user_defined_window(2)*100000)/100000;
R_squared=round((correlation_user_defined_window(2)*correlation_user_defined_window(2))*1000)/1000;
R_value_for_figure=strcat('Correlation Coefficient, R : ', mat2str(R_value));
R_squared_for_figure=strcat('Coefficient of Determination, R^2 : ', mat2str(R_squared));

%Save polyfit results to file to write out
Polyfit_summary_user(replicates1,1)=polyfit_x_value;
Polyfit_summary_user(replicates1,2)=polyfit_consent;
R_values_summary_user(replicates1,1)=R_value;
R_values_summary_user(replicates1,2)=R_squared;

%%Images of data within user defined window
%3D scatter plot showing the delta difference between replicates
Combined_delta2=(abs(Delta_height2))+(abs(Delta_width2))+abs(Delta_center2); %measure of total divergence between replicates
%Sorted Combined delta for use in figures
Combined_delta3=sort(Combined_delta2);

%Figure 3D scatter plot
image_3D_scatter_pr=figure;
image=scatter3(abs(Delta_center2), abs(Delta_height2), abs(Delta_width2), 50, log10(Combined_delta2),'filled');
colourbar_figure1= colorbar('location', 'SouthOutside');
xlabel(colourbar_figure1, 'Divergence');
% Add title and axis labels
title('Properties of Gaussian between replicates');
xlabel('Delta Center');
ylabel('Delta Height');
zlabel('Delta Width');
view(-45, 20);
threeD_scatter_plot_user_defined_window=strcat('Comparsion of the guassian properties between all proteins with one guassian between replicates',...
    mat2str(replicate_to_align_against),'_compared to_replicate',mat2str(replicates2),'_data_within_user_settings');

%generate image name
Image_name=strcat('image_',mat2str(Image_couter),'.png');

%Save image
saveas(image_3D_scatter_pr, Image_name);

%Save file name to image record
fprintf (Image_record,'%s,%s,\n',...
Image_name, threeD_scatter_plot_user_defined_window);

%Move output to Figure folder
movefile(Image_name, 'Alignment Figures');

%Add one to image counter
Image_couter=Image_couter+1;

%Figure line graphe
Alignment_graphic2= figure();
scatter(x_values4, y_values4, 50,  abs(z_values4), 'filled');
xlim([0,55]);
ylim([0,55]);
lsline
title({'Comparsion of proteins with single guassian between replicates'});
xlabel_name= strcat('Biological replicate',mat2str(replicate_to_align_against));
ylabel_name= strcat('Biological replicate',mat2str(replicates2));
xlabel(xlabel_name);
ylabel(ylabel_name);
%Add text with R, R2 and equation of the line
text(7, 48, Equation_of_polyfit);
text(7, 45, R_value_for_figure);
text(7, 42, R_squared_for_figure);
%Add a colorbar with tick labels
colourbar_figure= colorbar('location', 'EastOutside');
ylabel(colourbar_figure, 'Delta Center');
line_plot_user_defined_window=strcat('Comparsion of proteins with single guassian between replicates',mat2str(replicate_to_align_against),'_compared to_replicate',mat2str(replicates2),'_data_within_user_settings');

%generate image name
Image_name=strcat('image_',mat2str(Image_couter),'.png');

%Save image
saveas(Alignment_graphic2, Image_name);

%Save file name to image record
fprintf (Image_record,'%s,%s,\n',...
Image_name, line_plot_user_defined_window);

%Move output to Figure folder
movefile(Image_name, 'Alignment Figures');

%Add one to image counter
Image_couter=Image_couter+1;

%Hist of Center, height and width widthin user defined setting
fh2 = figure; % open figure tablet to create figure on
Ninty_precentile=round((Linear_list_lenght_user)*0.90); % determines the protein which is the top 90% precent

%hist for height
P1=subplot(4,1,1); %Define
number_of_precentile_height=Delta_height3(Ninty_precentile);
raw_bin_height_values_user=[0.1:0.1:7.5]; % bin values can be changed by the user to improve figure quility
Histc_height=histc(Delta_height3,raw_bin_height_values_user);
bar(raw_bin_height_values_user,Histc_height); %crete bar graph of values
hold on %keep bar graph and overlay plot
x_height=[number_of_precentile_height,number_of_precentile_height];
ylim_vlaues_P1=ylim(P1); %defined y limits from the bar graph
Max_height=ylim_vlaues_P1(2); %defined y limits max
Min_height=ylim_vlaues_P1(1); %defined y limits min
y_height=[Min_height,Max_height];
plot(x_height,y_height);
title(P1,'Height between replicates');
xlabel(P1,'Height');
ylabel(P1,'# of Gaussians');

%hist for center
P2=subplot(4,1,2); 
number_of_precentile_center=Delta_center3(Ninty_precentile);
raw_bin_center_values_user=[0.1:0.1:4.5];
Histc_center=histc(Delta_center3,raw_bin_center_values_user);
bar(raw_bin_center_values_user,Histc_center);
hold on %keep bar graph and overlay plot
x_center=[number_of_precentile_center,number_of_precentile_center];
ylim_vlaues_P2=ylim(P2); %defined y limits from the bar graph
Max_center=ylim_vlaues_P2(2); %defined y limits max
Min_center=ylim_vlaues_P2(1); %defined y limits min
y_center=[Min_center,Max_center];
plot(x_center,y_center);
title(P2,'Center between replicates');
xlabel(P2,'Center');
ylabel(P2,'# of Gaussians');

%hist for width
P3=subplot(4,1,3); 
number_of_precentile_width=Delta_width3(Ninty_precentile);
raw_bin_width_values_user=[0.1:0.1:4.5];
Histc_width=histc(Delta_width3,raw_bin_width_values_user);
bar(raw_bin_width_values_user,Histc_width);
hold on %keep bar graph and overlay plot
x_width=[number_of_precentile_width,number_of_precentile_width];
ylim_vlaues_P3=ylim(P3); %defined y limits from the bar graph
Max_width=ylim_vlaues_P3(2); %defined y limits max
Min_width=ylim_vlaues_P3(1); %defined y limits min
y_width=[Min_width,Max_width];
plot(x_width,y_width);
title(P3,'Width between replicates');
xlabel(P3,'Width');
ylabel(P3,'# of Gaussians');

%hist for Euclidean distance
P4=subplot(4,1,4); 
number_of_precentile_Euclidean_distance=Euclidean_distance3(Ninty_precentile);
Euclidean_distance_maximum_value=max(Euclidean_distance3);
raw_bin_Euclidean_distance_values_user=[0.5:0.5:Euclidean_distance_maximum_value];
Histc_Euclidean_distance=histc(Combined_delta3,raw_bin_Euclidean_distance_values_user);
bar(raw_bin_Euclidean_distance_values_user,Histc_Euclidean_distance);
hold on %keep bar graph and overlay plot
x_Euclidean_distance=[number_of_precentile_Euclidean_distance,number_of_precentile_Euclidean_distance];
ylim_vlaues_P4=ylim(P4); %defined y limits from the bar graph
Max_Euclidean_distance=ylim_vlaues_P4(2); %defined y limits max
Min_Euclidean_distance=ylim_vlaues_P4(1); %defined y limits min
y_Euclidean_distance=[Min_Euclidean_distance,Max_Euclidean_distance];
plot(x_Euclidean_distance,y_Euclidean_distance);
title(P4,'Euclidean Distance between replicates');
xlabel(P4,'Euclidean Distance');
ylabel(P4,'# of Gaussians');

%hist for combined, code to make hist of the divergents
%P4=subplot(4,1,4); 
%number_of_precentile_Combined=Combined_delta3(Ninty_precentile);
%raw_bin_combined_values_user=[0.1:0.25:12];
%Histc_combined=histc(Combined_delta3,raw_bin_combined_values_user);
%bar(raw_bin_combined_values_user,Histc_combined);
%hold on %keep bar graph and overlay plot
%x_combined=[number_of_precentile_Combined,number_of_precentile_Combined];
%ylim_vlaues_P4=ylim(P4); %defined y limits from the bar graph
%Max_combined=ylim_vlaues_P4(2); %defined y limits max
%Min_combined=ylim_vlaues_P4(1); %defined y limits min
%y_combined=[Min_combined,Max_combined];
%plot(x_combined,y_combined);
%title(P4,'Divergents between replicates');
%xlabel(P4,'Total Divergents');
%ylabel(P4,'# of Gaussians');

histo_plot_name_user_defined=strcat('Histogram of divergence between replicate',mat2str(replicate_to_align_against),'_compared to_replicate',mat2str(replicates2),'_data_within_user_settings');

%generate image name
Image_name=strcat('image_',mat2str(Image_couter),'.png');

%Save image
saveas(fh2, Image_name);
close(fh2);

%Save file name to image record
fprintf (Image_record,'%s,%s,\n',...
Image_name, histo_plot_name_user_defined);

%Move output to Figure folder
movefile(Image_name, 'Alignment Figures');

%Add one to image counter
Image_couter=Image_couter+1;


end
end

%%Alignment output
%Create a table with polynomial fit adjustmetns for replicates
%corresponding to datapoint within the user define value
Summary_poly_fit=zeros(Unique_replicate,2);
postions_to_write=length(replicates_minus_alignment);
Summary_poly_fit(y_maximum,1)=1;
Summary_poly_fit(y_maximum,2)=0;
Summary_R_values(y_maximum,1)=1;
Summary_R_values(y_maximum,2)=1;
for counter_summary= 1:postions_to_write
    Summary_poly_fit(replicates_minus_alignment(counter_summary),1)= Polyfit_summary_user(counter_summary,1);
    Summary_poly_fit(replicates_minus_alignment(counter_summary),2)= Polyfit_summary_user(counter_summary,2);
    Summary_R_values(replicates_minus_alignment(counter_summary),1)= R_values_summary_user(counter_summary,1);
    Summary_R_values(replicates_minus_alignment(counter_summary),2)= R_values_summary_user(counter_summary,2);
end

%Write out table of polynomial fit adjustmetns for replicates
fid6A = fopen('Summary_table_polynomial_fit_adjustmetns_for_replicates_user_defined_thresholds.csv','w');
fprintf (fid6A,'%s,%s,%s,%s,%s\n',...
'Replicate Number', 'Indeterminates', 'Coefficients', 'R value', 'R squared');  %Write Header
for ix=1:Unique_replicate
    fprintf(fid6A, '%6.4f,%6.4f,%6.4f,%6.4f,%6.4f\n',...
        (ix), Summary_poly_fit(ix, 1), Summary_poly_fit(ix, 2), Summary_R_values(ix,1),Summary_R_values(ix,2));
end
fclose(fid6A);


%Create a table with polynomial fit adjustmetns for replicates
%corresponding to datapoint within the user define value
Summary_poly_fit1=zeros(Unique_replicate,2);
postions_to_write=length(replicates_minus_alignment);
Summary_poly_fit1(y_maximum,1)=1;
Summary_poly_fit1(y_maximum,2)=0;
Summary_R_values1(y_maximum,1)=1;
Summary_R_values1(y_maximum,2)=1;
for counter_summary= 1:postions_to_write
    Summary_poly_fit1(replicates_minus_alignment(counter_summary),1)= Polyfit_summary_global(counter_summary,1);
    Summary_poly_fit1(replicates_minus_alignment(counter_summary),2)= Polyfit_summary_global(counter_summary,2);
    Summary_R_values1(replicates_minus_alignment(counter_summary),1)= R_values_summary_global(counter_summary,1);
    Summary_R_values1(replicates_minus_alignment(counter_summary),2)= R_values_summary_global(counter_summary,2);
end

%Write out table of polynomial fit adjustmetns for replicates
fid6B = fopen('Summary_table_polynomial_fit_adjustmetns_for_replicates_global.csv','w');
fprintf (fid6B,'%s,%s,%s,%s,%s\n',...
'Replicate Number', 'Indeterminates', 'Coefficients', 'R value', 'R squared');  %Write Header
for ix=1:Unique_replicate
    fprintf(fid6B, '%6.4f,%6.4f,%6.4f,%6.4f,%6.4f\n',...
        (ix), Summary_poly_fit1(ix, 1), Summary_poly_fit1(ix, 2), Summary_R_values1(ix,1),Summary_R_values1(ix,2));
end
fclose(fid6B);

%Move output to Realignment folder
movefile('Summary_table_polynomial_fit_adjustmetns_for_replicates_user_defined_thresholds.csv', 'Realignment');
movefile('Summary_table_polynomial_fit_adjustmetns_for_replicates_global.csv', 'Realignment');
     
%%%Using Alignment output to correct gaussians 

%%Using the alignment data to recalibrate replicates
%adjust center of Combined_OutputGaus by equation determined by polyfit
%Note: currently designed to use the MvsL as is generated first, if the
%user requires HvsL modify the experimental channel loop
if isequal(Alignment_to_user, Experimental_channel)
cd('Realignment');
Alignment_correction_values=importdata('Summary_table_polynomial_fit_adjustmetns_for_replicates_user_defined_thresholds.csv');
cd('../');
end

%% Using fitted curve adjust observed replicate data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%y_maximum is the replicate which the other replicates are aligned to
replicate_to_align_against= y_maximum;
    
%impott list of proteins with a single guassian from replicates to test
Gaus_import_Name_for_adjustment_reference= strcat('Processed Gaussian\',Experimental_channel,'_Combined_OutputGaus_rep',mat2str(replicate_to_align_against),'.csv');
Gaus_import_reference_rep=importdata(Gaus_import_Name_for_adjustment_reference);
Protein_names_gaussian_for_reference_length=length(Gaus_import_reference_rep.data);
Protein_names_gaussian_for_reference=Gaus_import_reference_rep.textdata(2:Protein_names_gaussian_for_reference_length+1,4);


%Copy raw data across to use for alignment
Gaus_import_reference_rep_raw=Gaus_import_reference_rep;

for alignment_counter=1:Unique_replicate
    alignment_indeterminates=Alignment_correction_values.data(alignment_counter,2);
    alignment_coefficients=Alignment_correction_values.data(alignment_counter,3);
    
%Import Matlab processed data
%Data to be adjusted based on polyfit equation
Gaus_import_Name_for_adjustment= strcat('Processed Gaussian\',Experimental_channel,'_Combined_OutputGaus_rep',mat2str(alignment_counter),'.csv');
Gaus_import_for_adjustment=importdata(Gaus_import_Name_for_adjustment);           %Import the Gaus data from previous scripts
Protein_names_gaussian_for_adjustment_length=length(Gaus_import_for_adjustment.data);
Protein_names_gaussian_for_adjustment=Gaus_import_for_adjustment.textdata(2:Protein_names_gaussian_for_adjustment_length+1,4);

Final_Gaus_import_for_adjustment=Gaus_import_for_adjustment;

%adjust raw data based on alignment
Raw_data_Name_for_adjustment= strcat('Processed Gaussian\',Experimental_channel,'_Raw_data_maxquant_rep',mat2str(alignment_counter),'.csv');
Raw_data_import_for_adjustment=importdata(Raw_data_Name_for_adjustment);           %Import the Gaus data from previous scripts
Raw_data_for_adjustment_size=size(Raw_data_import_for_adjustment.data);

%adjust raw data based on alignment
HvsM_Raw_data_Name_for_adjustment= strcat('Processed Gaussian\HvsM_Raw_data_maxquant_rep',mat2str(alignment_counter),'.csv');
HvsM_Raw_data_import_for_adjustment=importdata(HvsM_Raw_data_Name_for_adjustment);           %Import the Gaus data from previous scripts
HvsM_Raw_data_for_adjustment_size=size(HvsM_Raw_data_import_for_adjustment.data);

%Copy raw data across to use for alignment
Gaus_import_for_adjustment_raw=Gaus_import_for_adjustment;
    
%Create array for storing gaussian information
Global_aligned_gaus.Protein_name_raw=cell(0,0);
Global_aligned_gaus.Center_replicate1_raw=zeros(0,0);
Global_aligned_gaus.Center_replicate2_being_tested_raw=zeros(0,0);
Global_aligned_gaus.Delta_center_raw=zeros(0,0);    
    
%Find Gaussians fits common to both replicates
for protein_counter=1:Protein_names_gaussian_for_adjustment_length
%find position of proteins
[replicate_location]=ind2sub(Protein_names_gaussian_for_reference_length,...
strmatch(Protein_names_gaussian_for_adjustment{protein_counter}, Protein_names_gaussian_for_reference(:),  'exact'));

%Define varibles to populate
number_guassian=length(replicate_location);
Center_to_test=zeros(number_guassian,1);

for protein_counter2=1:number_guassian 
    
%find the corresponding center
Center_to_test(protein_counter2,1)=Gaus_import_reference_rep_raw.data(replicate_location(protein_counter2),2);

end

if ~isempty(Center_to_test)
    
%find position of closest value to expected center
[row_1 row_1]=min(abs(Center_to_test- Gaus_import_for_adjustment_raw.data(protein_counter,2)));

%test if the cloest value is within user defined window
if abs(Gaus_import_for_adjustment_raw.data(protein_counter,2) - Gaus_import_reference_rep_raw.data(replicate_location(row_1),2)) < User_alignment_window1
    
Global_aligned_gaus.Protein_name_raw(length(Global_aligned_gaus.Protein_name_raw)+1)...
=Protein_names_gaussian_for_adjustment(protein_counter);
Global_aligned_gaus.Center_replicate1_raw(length(Global_aligned_gaus.Center_replicate1_raw)+1)...
= Gaus_import_for_adjustment_raw.data(protein_counter,2);
Global_aligned_gaus.Center_replicate2_being_tested_raw(length(Global_aligned_gaus.Center_replicate2_being_tested_raw)+1)...
= Gaus_import_reference_rep_raw.data(replicate_location(row_1),2);
Global_aligned_gaus.Delta_center_raw(length(Global_aligned_gaus.Delta_center_raw)+1)...
=abs(Gaus_import_for_adjustment_raw.data(protein_counter,2) - Gaus_import_reference_rep_raw.data(replicate_location(row_1),2));

end
end
end    
    
    
%% alignment within user defined values
%Create array for storing gaussian information
Global_aligned_gaus_1st_fit.Protein_name=cell(0,0);
Global_aligned_gaus_1st_fit.Center_replicate1=zeros(0,0);
Global_aligned_gaus_1st_fit.Center_replicate2_being_tested=zeros(0,0);
Global_aligned_gaus_1st_fit.Delta_center=zeros(0,0);
   
%generate new center values for adjusted replicate using just coefficent
Gaus_import_for_adjustment_1st=Gaus_import_for_adjustment;
Gaus_import_for_adjustment_1st.data(:,2)=(Gaus_import_for_adjustment_1st.data(:,2)-alignment_coefficients)/alignment_indeterminates;
%Adjust height accordingly
Gaus_import_for_adjustment_1st.data(:,1)=(Gaus_import_for_adjustment_1st.data(:,1)*alignment_indeterminates);
%Adjust weight accordingly
Gaus_import_for_adjustment_1st.data(:,3)=(Gaus_import_for_adjustment_1st.data(:,3)*1/alignment_indeterminates);
Gaus_import_for_adjustment_1st.data(:,6)=(round((SEC_fit(1)*Gaus_import_for_adjustment_1st.data(:,2)+SEC_fit(2))/10)*10);

%Find Gaussians fits common to both replicates
for protein_counter=1:Protein_names_gaussian_for_adjustment_length
%find position of proteins
[replicate_location]=ind2sub(Protein_names_gaussian_for_reference_length,...
strmatch(Protein_names_gaussian_for_adjustment{protein_counter}, Protein_names_gaussian_for_reference(:),  'exact'));

%Define varibles to populate
number_guassian=length(replicate_location);
Center_to_test=zeros(number_guassian,1);

for protein_counter2=1:number_guassian 
    
%find the corresponding center
Center_to_test(protein_counter2,1)=Gaus_import_reference_rep.data(replicate_location(protein_counter2),2);

end

if ~isempty(Center_to_test)
    
%find position of closest value to expected center
[row_1 row_1]=min(abs(Center_to_test- Gaus_import_for_adjustment_1st.data(protein_counter,2)));

%test if the cloest value is within user defined window
if abs(Gaus_import_for_adjustment_1st.data(protein_counter,2) - Gaus_import_reference_rep.data(replicate_location(row_1),2)) < User_alignment_window2
    
Global_aligned_gaus_1st_fit.Protein_name(length(Global_aligned_gaus_1st_fit.Protein_name)+1)...
=Protein_names_gaussian_for_adjustment(protein_counter);
Global_aligned_gaus_1st_fit.Center_replicate1(length(Global_aligned_gaus_1st_fit.Center_replicate1)+1)...
= Gaus_import_for_adjustment_1st.data(protein_counter,2);
Global_aligned_gaus_1st_fit.Center_replicate2_being_tested(length(Global_aligned_gaus_1st_fit.Center_replicate2_being_tested)+1)...
= Gaus_import_reference_rep.data(replicate_location(row_1),2);
Global_aligned_gaus_1st_fit.Delta_center(length(Global_aligned_gaus_1st_fit.Delta_center)+1)...
=abs(Gaus_import_for_adjustment_1st.data(protein_counter,2) - Gaus_import_reference_rep.data(replicate_location(row_1),2));

end
end
end


%% Using crude alined data fit based on density

xi=linspace(min(Global_aligned_gaus_1st_fit.Center_replicate1(:)),max(Global_aligned_gaus_1st_fit.Center_replicate1(:)),fraction_number(2));
yi=linspace(min(Global_aligned_gaus_1st_fit.Center_replicate2_being_tested(:)),max(Global_aligned_gaus_1st_fit.Center_replicate2_being_tested(:)),fraction_number(2));

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

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

%find position max with within each column
Values_to_correct=zeros(0,0);

%create array 
Maximum_position=zeros(1,(fraction_number(2)));

for maximum_counter=1:(fraction_number(2))
%find the maximum value
idx = max(Z(maximum_counter,:));
if ~idx==0 && ~isempty(idx)==1
position=find(Z(:,maximum_counter),idx);
number_of_positions=length(position);
sum_of_position=sum(position);
Maximum_position(maximum_counter)=sum_of_position/number_of_positions;
elseif idx==0 | isempty(idx)==1
    Values_to_correct(length(Values_to_correct)+1,1)=maximum_counter;
    Maximum_position(maximum_counter)=0;    
end
end

%Correct any NaN values   
Maximum_position(isnan(Maximum_position))=0;

%Fix missing values
for fixing_value_counter=1:length(Values_to_correct)
    
    %find number before missing value with measurement
    ready1=false;
    count_backwards=1;

    while ready1==false
        if ~Maximum_position(Values_to_correct(fixing_value_counter)-count_backwards) ==0
        position_1_value=Maximum_position(Values_to_correct(fixing_value_counter)-count_backwards);
        ready1=true;
        else
        count_backwards=count_backwards+1;
        end
    end
    
    %find number after missing value with measurement
    ready2=false;
    count_forwards=1;
    
    while ready2==false
        if ~Maximum_position(Values_to_correct(fixing_value_counter)+count_forwards) ==0
        position_2_value=Maximum_position(Values_to_correct(fixing_value_counter)+count_forwards);
        ready2=true;
        else
        count_forwards=count_forwards+1;
        end
    end
    
    %Fit line to values to determine what the missing values is
    add_to_missing_value=round((position_2_value-position_1_value)/(count_forwards+count_backwards));
  
    %Write value to Maximum position
    Maximum_position(Values_to_correct(fixing_value_counter))=position_1_value+add_to_missing_value;
end
    

%Create alignment based on density
smoothed_maximum_postions=smooth(Maximum_position,3);

%Create alignment template 
aligned_positions=1:1:(fraction_number(2))';
Corrected_positions= smoothed_maximum_postions;

%Determine equation of the density line and correlation coeffectinces
Density_polyfit=polyfit(aligned_positions',Corrected_positions,1);
Density_correlation=corrcoef(aligned_positions',Corrected_positions);

%Save values to be written out
Density_polyfit_to_write_out(alignment_counter,1)=Density_polyfit(1);
Density_polyfit_to_write_out(alignment_counter,2)=Density_polyfit(2);
Density_correlation_to_write_out(alignment_counter,1)=round(Density_correlation(2)*100000)/100000;
Density_correlation_to_write_out(alignment_counter,2)=round((Density_correlation(2)*Density_correlation(2))*1000)/1000;

%for visually checking alignment
%surf(Z) % Use Surf to look at density plot of alignment
%plot(1:fraction_number(2),smooth(Maximum_position,3)); % Use plot to look at density plot of alignment
%Create Figure of all the alignment data


%% alignment within user defined values
%Create array for storing gaussian information
Global_aligned_gaus_2ndfit.Protein_name=cell(0,0);
Global_aligned_gaus_2ndfit.Center_replicate1=zeros(0,0);
Global_aligned_gaus_2ndfit.Center_replicate2_being_tested=zeros(0,0);
Global_aligned_gaus_2ndfit.Delta_center=zeros(0,0);

%Combined alignment based on density and single protein gaussians
Final_realignment_value(alignment_counter,1)=Density_polyfit_to_write_out(alignment_counter,1)*Alignment_correction_values.data(alignment_counter,2);
Final_realignment_value(alignment_counter,2)=Density_polyfit_to_write_out(alignment_counter,2)+Alignment_correction_values.data(alignment_counter,3);

%generate new center values for adjusted replicate using just coefficent
Gaus_import_for_adjustment_2nd=Gaus_import_for_adjustment_1st;
Gaus_import_for_adjustment_2nd.data(:,2)=(Gaus_import_for_adjustment_1st.data(:,2)-Final_realignment_value(alignment_counter,2))/Final_realignment_value(alignment_counter,1);
Gaus_import_for_adjustment_2nd.data(:,6)=(round((SEC_fit(1)*Gaus_import_for_adjustment_2nd.data(:,2)+SEC_fit(2))/10)*10);
   
%Find Gaussians fits common to both replicates
for protein_counter=1:Protein_names_gaussian_for_adjustment_length
%find position of proteins
[replicate_location]=ind2sub(Protein_names_gaussian_for_reference_length,...
strmatch(Protein_names_gaussian_for_adjustment{protein_counter}, Protein_names_gaussian_for_reference(:),  'exact'));

%Define varibles to populate
number_guassian=length(replicate_location);
Center_to_test=zeros(number_guassian,1);

for protein_counter2=1:number_guassian 
    
%find the corresponding center
Center_to_test(protein_counter2,1)=Gaus_import_reference_rep.data(replicate_location(protein_counter2),2);

end

if ~isempty(Center_to_test)
    
%find position of closest value to expected center
[row_1 row_1]=min(abs(Center_to_test- Gaus_import_for_adjustment.data(protein_counter,2)));

%test if the cloest value is within dyanamic user defined window
if abs(Gaus_import_for_adjustment.data(protein_counter,2) - Final_realignment_value(alignment_counter,1)*Gaus_import_reference_rep.data(replicate_location(row_1),2)+Final_realignment_value(alignment_counter,2)) < (User_alignment_window2) 
% by using a dynamic window for the final realignment replicates with a different slopes can still be compared
    
Global_aligned_gaus_2ndfit.Protein_name(length(Global_aligned_gaus_2ndfit.Protein_name)+1)...
=Protein_names_gaussian_for_adjustment(protein_counter);
Global_aligned_gaus_2ndfit.Center_replicate1(length(Global_aligned_gaus_2ndfit.Center_replicate1)+1)...
= Gaus_import_for_adjustment.data(protein_counter,2);
Global_aligned_gaus_2ndfit.Center_replicate2_being_tested(length(Global_aligned_gaus_2ndfit.Center_replicate2_being_tested)+1)...
= Gaus_import_reference_rep.data(replicate_location(row_1),2);
Global_aligned_gaus_2ndfit.Delta_center(length(Global_aligned_gaus_2ndfit.Delta_center)+1)...
=abs(Gaus_import_for_adjustment.data(protein_counter,2) - Gaus_import_reference_rep.data(replicate_location(row_1),2));

end
end
end

Final_polyfit=polyfit(Global_aligned_gaus_2ndfit.Center_replicate2_being_tested,Global_aligned_gaus_2ndfit.Center_replicate1,1);
Final_correlation=corrcoef(Global_aligned_gaus_2ndfit.Center_replicate1,Global_aligned_gaus_2ndfit.Center_replicate2_being_tested);

%Save values to be written out
Final_polyfit_to_write_out(alignment_counter,1)=Final_polyfit(1);
Final_polyfit_to_write_out(alignment_counter,2)=Final_polyfit(2);
Final_correlation_to_write_out(alignment_counter,1)=round(Final_correlation(2)*100000)/100000;
Final_correlation_to_write_out(alignment_counter,2)=round((Final_correlation(2)*Final_correlation(2))*1000)/1000;

%Write out table of polynomial fit tables
cd('Realignment');
    fprintf(fid6C, '%6.4f,%6.4f,%6.4f,%6.4f,%6.4f\n',...
        (alignment_counter), Density_polyfit_to_write_out(alignment_counter, 1), Density_polyfit_to_write_out(alignment_counter, 2),...
        Density_correlation_to_write_out(alignment_counter,1),Density_correlation_to_write_out(alignment_counter,2));

    fprintf(fid6D, '%6.4f,%6.4f,%6.4f,%6.4f,%6.4f\n',...
        (alignment_counter), Final_polyfit_to_write_out(alignment_counter, 1), Final_polyfit_to_write_out(alignment_counter, 2),...
        Final_correlation_to_write_out(alignment_counter,1),Final_correlation_to_write_out(alignment_counter,2));
cd('../');

%Calulate what each fraction corresponds to with the new scale
%%Using the alignment data to recalibrate replicates
%adjust center of Combined_OutputGaus by equation determined by polyfit
%Note: currently designed to use the MvsL as is generated first, if the
%user requires HvsL modify the experimental channel loop
if isequal(Alignment_to_user, Experimental_channel)
cd('Realignment');
FinaL_alignment_correction=importdata('Summary_table_polynomial_fit_both_adjustmetns_for_replicates.csv');
cd('../');
end

%Generate array with fractions -5 and +60 to enable easy realignment of experiments
New_Fraction_numbering=[-5:60];
lenght_new_fraction=length(New_Fraction_numbering);

%Create array to write to
adjusted_raw_data=zeros(Raw_data_for_adjustment_size(1),lenght_new_fraction);
HvsM_adjusted_raw_data=zeros(Raw_data_for_adjustment_size(1),lenght_new_fraction);

% linear calibration step based the results of the linear and density
% alignments

re_calibrated_fraction_number=([1:fraction_number(2)]-FinaL_alignment_correction.data(alignment_counter,3))/FinaL_alignment_correction.data(alignment_counter,2);
maximum_fraction_number=round(max(re_calibrated_fraction_number));

%only re align files that need to be adjusted
if  ~(FinaL_alignment_correction.data(alignment_counter,2)==1 && FinaL_alignment_correction.data(alignment_counter,3)==0)

%determine the first fraction to write to
First_position = abs(New_Fraction_numbering-re_calibrated_fraction_number(1));
[First_position_values First_position_values]= min(First_position);  


for raw_data_adjustment1=1:Raw_data_for_adjustment_size(1)
    
    %Reset the writing position back to first fraction  
    writing_position=First_position_values;
    
    for raw_data_adjustment2= 2:maximum_fraction_number
           
    %find position in the array flanking the fraction to adjust to 
    position_1 = max(find(re_calibrated_fraction_number(:)<(raw_data_adjustment2-1)));
    position_2 = min(find(re_calibrated_fraction_number(:)>(raw_data_adjustment2-1)));
    
    %Correction value
    position_1_values=abs(abs(re_calibrated_fraction_number(position_1))-ceil(re_calibrated_fraction_number(position_1)));
    position_2_values=1-abs(abs(re_calibrated_fraction_number(position_2))-ceil(re_calibrated_fraction_number(position_2)));
   
    
    %Check position values contains values, if it does demostrates
    %fractions missed at the start or end of the run
    %check position 1
    if isempty(position_1_values)
        position_1_values=0;
    end
    
    %check position 2
    if isempty(position_2_values)
        position_2_values=0;
    end
    
    %determine the precentage contribution from each value
    sum_value_p1_p2=position_1_values+position_2_values;
    
    %Create value to write average to    
    raw_values_readjusted=zeros(1,2);
    HvsM_raw_values_readjusted=zeros(1,2);
    
    %test if position contains a isotopologue ratio
    if ~(position_1_values == 0 & position_2_values == 0)
    
        %Determine positions to sum
        if  isempty(position_1)
         
        raw_values_readjusted(1)=(1-position_2_values)*Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_2+1); % plus 1 for position as array contains replicate number
        raw_values_readjusted(2)=position_2_values*Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_2+1); % plus 1 for position as array contains replicate number

        %HvsM adjustment
        HvsM_raw_values_readjusted(1)=(1-position_2_values)*HvsM_Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_2+1);  % plus 1 for position as array contains replicate number
        HvsM_raw_values_readjusted(2)=position_2_values*HvsM_Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_2+1);  % plus 1 for position as array contains replicate number
        
        elseif ~(isempty(position_1) & isempty(position_2))
 
        raw_values_readjusted(1)=(position_1_values/sum_value_p1_p2)*Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_1+1);      
        raw_values_readjusted(2)=(position_2_values/sum_value_p1_p2)*Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_2+1);

        %HvsM adjustment
        HvsM_raw_values_readjusted(1)=(position_1_values/sum_value_p1_p2)*HvsM_Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_1+1);
        HvsM_raw_values_readjusted(2)=(position_2_values/sum_value_p1_p2)*HvsM_Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_2+1);

        elseif isempty(position_2)
            
        raw_values_readjusted(1)=position_1_values*Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_1+1); % plus 1 for position as array contains replicate number
        raw_values_readjusted(2)=(1-position_1_values)*Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_1+1); % plus 1 for position as array contains replicate number

        %HvsM adjustment
        HvsM_raw_values_readjusted(1)=position_1_values*HvsM_Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_1+1);  % plus 1 for position as array contains replicate number
        HvsM_raw_values_readjusted(2)=(1-position_1_values)*HvsM_Raw_data_import_for_adjustment.data(raw_data_adjustment1,position_1+1);  % plus 1 for position as array contains replicate number    

        end
     
    %remove any empty values
    NaN_check=isnan(raw_values_readjusted);
    HvsM_NaN_check=isnan(HvsM_raw_values_readjusted);
    raw_values_readjusted(NaN_check)=[];  
    HvsM_raw_values_readjusted(HvsM_NaN_check)=[]; 

    %Add the adjusted value to the raw_data_array
    adjusted_raw_data(raw_data_adjustment1,writing_position)= sum(raw_values_readjusted);
    HvsM_adjusted_raw_data(raw_data_adjustment1,writing_position)= sum(HvsM_raw_values_readjusted);
    
    elseif (position_1_values == 0 & position_2_values == 0)
        
    adjusted_raw_data(raw_data_adjustment1,writing_position)=0;
    HvsM_adjusted_raw_data(raw_data_adjustment1,writing_position)=0;
    
    end
    %counter to write to
    writing_position=writing_position+1;
    
    end   
end
  

elseif (FinaL_alignment_correction.data(alignment_counter,2)==1 && FinaL_alignment_correction.data(alignment_counter,3)==0)
        
        %determine the first fraction to write to
        First_position = abs(New_Fraction_numbering-re_calibrated_fraction_number(1));
        [First_position_values First_position_values]= min(First_position); 

        %Copy data not needing to be adjusted directly into adjusted_raw_data array
        for raw_data_adjustment1=1:Raw_data_for_adjustment_size(1) 
                  
        writing_position=First_position_values;

        for raw_data_adjustment2=2:Raw_data_for_adjustment_size(2)
    
        adjusted_raw_data(raw_data_adjustment1,writing_position)=Raw_data_import_for_adjustment.data(raw_data_adjustment1,raw_data_adjustment2);
        HvsM_adjusted_raw_data(raw_data_adjustment1,writing_position)=HvsM_Raw_data_import_for_adjustment.data(raw_data_adjustment1,raw_data_adjustment2);

        writing_position=writing_position+1;
        
        end
        end

end

%Remove zeros and replace with NaN
adjusted_raw_data(adjusted_raw_data==0)=NaN;
HvsM_adjusted_raw_data(HvsM_adjusted_raw_data==0)=NaN;

%Readjust Gaussians based on Final alignment
Gaus_import_for_adjustment.data(:,2)=(Final_Gaus_import_for_adjustment.data(:,2)-FinaL_alignment_correction.data(alignment_counter,3))/FinaL_alignment_correction.data(alignment_counter,2);
%Adjust height accordingly
Gaus_import_for_adjustment.data(:,1)=(Gaus_import_for_adjustment.data(:,1)*FinaL_alignment_correction.data(alignment_counter,2));
%Adjust weight accordingly
Gaus_import_for_adjustment.data(:,3)=(Gaus_import_for_adjustment.data(:,3)*1/FinaL_alignment_correction.data(alignment_counter,2));
Gaus_import_for_adjustment.data(:,6)=(round((SEC_fit(1)*Final_Gaus_import_for_adjustment.data(:,2)+SEC_fit(2))/10)*10);


%adjusted replicate chromotgram
Adjusted_Chromatograms=zeros(Protein_names_gaussian_for_adjustment_length,Chromatograms_axis);
for Counter_eleven= 1:Protein_names_gaussian_for_adjustment_length
    value_for_adjusted_Height=Gaus_import_for_adjustment.data(Counter_eleven,1);
    value_for_adjusted_Center=Gaus_import_for_adjustment.data(Counter_eleven,2);
    value_for_adjusted_Width=Gaus_import_for_adjustment.data(Counter_eleven,3);    
    for Counter_twelve =1:Chromatograms_axis
    Adjusted_Chromatograms(Counter_eleven,Counter_twelve)...
        =  value_for_adjusted_Height*exp(-((Counter_twelve-(value_for_adjusted_Center+5))/value_for_adjusted_Width).^2);
    end
end


%Generate Combined Gaussian curves for euclidean distisnce calculations
Combining_Adjusted_Chromatograms=zeros(Protein_names_gaussian_for_adjustment_length,Chromatograms_axis);
for combing_counter=1:Protein_names_gaussian_for_adjustment_length
    find_string_combining= (strmatch(Protein_names_gaussian_for_adjustment{combing_counter},Gaus_import_for_adjustment.textdata(:,4), 'exact')-1);
    number_of_string_to_combine=length(find_string_combining);
    for internal_combine_counter=1:number_of_string_to_combine
        Combining_Adjusted_Chromatograms(combing_counter,:)= Combining_Adjusted_Chromatograms(combing_counter,:) + Adjusted_Chromatograms(find_string_combining(internal_combine_counter),:);
    end
end

%Remove near zero values with a user defined value
User_defined_zero_value=0.01; %lowest value to be shown in Adjusted Chromatograms
Combining_Adjusted_Chromatograms(Combining_Adjusted_Chromatograms<0.1) = User_defined_zero_value(1);
%The user can adjust the lowest value reported in the adjusted
%chromatogram, suggested 0.1

%Write out adjusted figures
%Figure of correlation between centers of raw data
Alignment_figure=figure;
Alignment_figure1=subplot(3,1,1);
%panel one, plot all raw data with linear fit
scatter(Global_aligned_gaus.Center_replicate2_being_tested_raw,...
Global_aligned_gaus.Center_replicate1_raw,4,...
abs(Global_aligned_gaus.Delta_center_raw), 'filled');
title({'Comparsion of single guassian fits between replicates'});
xlabel_name= strcat('Biological replicate',mat2str(alignment_counter));
ylabel_name= strcat('Biological replicate',mat2str(replicate_to_align_against));
xlabel(xlabel_name); 
ylabel(ylabel_name);
%plot fitted lined
x_plot_value=1:1:55;
y_plot_value=1:1:55;
corrected_x_plot_value=(x_plot_value-alignment_coefficients)/alignment_indeterminates;
hold on;
plot(corrected_x_plot_value,y_plot_value,'-','MarkerSize', 4, 'Color', [0.9 0.2 .4]);
%Add colourbar
colourbar_figure= colorbar('location', 'EastOutside');
ylabel(colourbar_figure, 'Delta Center');
%Add limits on figure
xlim_max=max(Global_aligned_gaus.Center_replicate2_being_tested_raw);
ylim_max=max(Global_aligned_gaus.Center_replicate1_raw);
xlim([0,round(xlim_max*1.02)]);
ylim([0,round(ylim_max*1.02)]);
Equation_of_line1=strcat('Equation of the fit:',mat2str(alignment_indeterminates),'x+',mat2str(alignment_coefficients));
%Add text
text(1, ylim_max*0.95, Equation_of_line1);
hold off; % reset hold on


%panel two, plot linear corrected values, within user defined window with
%density based fit  %%%%chnage to have densitycoloured
Alignment_figure2=subplot(3,1,2);
%panel one, plot all raw data with linear fit
scatter(Global_aligned_gaus_1st_fit.Center_replicate2_being_tested,...
    Global_aligned_gaus_1st_fit.Center_replicate1,8,...
abs(Global_aligned_gaus_1st_fit.Delta_center), 'filled');
title({'Comparsion of all guassian fits within 2 fractions between replicates'});
xlabel_name= strcat('Biological replicate',mat2str(alignment_counter));
ylabel_name= strcat('Biological replicate',mat2str(replicate_to_align_against));
xlabel(xlabel_name); 
ylabel(ylabel_name);
%Add colourbar
colourbar_figure= colorbar('location', 'EastOutside');
ylabel(colourbar_figure, 'Delta Center');
%Add limits on figure
hold on
%Add density based plot over values
corrected_x_plot_value1=(x_plot_value-Density_polyfit_to_write_out(alignment_counter,2))/Density_polyfit_to_write_out(alignment_counter,1);
hold on;
plot(corrected_x_plot_value1,y_plot_value,'-','MarkerSize', 4, 'Color', [0.9 0.2 .4]);
%Add limits on figure
xlim_max=max(Global_aligned_gaus_1st_fit.Center_replicate2_being_tested);
ylim_max=max(Global_aligned_gaus_1st_fit.Center_replicate1);
caxis_max=max(abs(Global_aligned_gaus_1st_fit.Delta_center));
xlim([0,round(xlim_max*1.02)]);
ylim([0,round(ylim_max*1.02)]);
caxis([0, round(caxis_max*1.1)]);
Equation_of_line2=strcat('Equation of the fit:',mat2str(round(Density_polyfit_to_write_out(alignment_counter,1)*1000)/1000),...
    'x+',mat2str(round(Density_polyfit_to_write_out(alignment_counter,2)*1000)/1000));
Equation_of_line2A=strcat('Number of Gaussians:',mat2str(length(Global_aligned_gaus_1st_fit.Center_replicate1)));
%Add text
text(1, ylim_max*0.95, Equation_of_line2);
text(1, ylim_max-15, Equation_of_line2A);
hold off % reset hold on

%panel two, plot linear corrected values, within user defined window with
%density based fit  %%%%chnage to have densitycoloured
Alignment_figure3=subplot(3,1,3);
%panel one, plot all raw data with linear fit
scatter(Global_aligned_gaus_2ndfit.Center_replicate2_being_tested,...
Global_aligned_gaus_2ndfit.Center_replicate1,8,...
abs(Global_aligned_gaus_2ndfit.Delta_center), 'filled');
title({'Final alignment of Gaussians'});
ylabel_name= strcat('Biological replicate',mat2str(replicate_to_align_against));
xlabel_name= strcat('Biological replicate',mat2str(alignment_counter));
xlabel(xlabel_name); 
ylabel(ylabel_name);
%Add colourbar
colourbar_figure= colorbar('location', 'EastOutside');
ylabel(colourbar_figure, 'Delta Center');
%Add limits on figure
xlim_max=max(Global_aligned_gaus_2ndfit.Center_replicate2_being_tested);
ylim_max=max(Global_aligned_gaus_2ndfit.Center_replicate1);
caxis_max=max(abs(Global_aligned_gaus_2ndfit.Delta_center));
xlim([0,round(xlim_max*1.02)]);
ylim([0,round(ylim_max*1.02)]);
caxis([0, round(caxis_max*1.1)])
Equation_of_line3=strcat('Combined fit:',mat2str(round(Final_polyfit_to_write_out(alignment_counter,1)*1000)/1000),...
    'x+',mat2str(round(Final_polyfit_to_write_out(alignment_counter,2)*1000)/1000));
Equation_of_line3A=strcat('Number of Gaussians:',mat2str(length(Global_aligned_gaus_2ndfit.Center_replicate1)));
%Add text
text(1, ylim_max*0.95, Equation_of_line3);
text(1, ylim_max-15, Equation_of_line3A);

%set size of subplots
s1Pos = [0.1300    0.7458    0.6512    0.1544];
s2Pos = [0.1300    0.4462    0.6512    0.1544];
s3Pos = [0.1300    0.1100    0.6512    0.2157];
set(Alignment_figure1,'position',s1Pos);
set(Alignment_figure2,'position',s2Pos);
set(Alignment_figure3,'position',s3Pos);

Alignment_name_to_record=strcat('Comparsion of alignment steps comparing replicate ',mat2str(replicate_to_align_against),' to replicate',mat2str(alignment_counter));

%generate image name
Image_name=strcat('image_',mat2str(Image_couter),'.png');
%Add one to image counter
Image_couter=Image_couter+1;

%Save file name to image record
fprintf (Image_record,'%s,%s,\n',...
Image_name, Alignment_name_to_record);

%Save image
cd('Alignment Figures');
saveas(Alignment_figure, Image_name);
cd('../')

%Write out adjusted chromotograms
cd('Realignment');

%Write out summary file 
fid_summary_1_name=strcat('fraction_overlap summary_',Experimental_channel,'_replicate_',mat2str(alignment_counter),'.csv');
fid_summary_1 = fopen(fid_summary_1_name,'at');
Experiment_name1=strcat('PCP_SEC_SILAC_',Experimental_channel,'_replicate_',mat2str(alignment_counter),'_uncalibrated fraction number');

%Write out header for summary file
Column_header=repmat('%s,', 1, (Raw_data_for_adjustment_size(2)));
for column_counter=1:Raw_data_for_adjustment_size(2)-1
    Column_name1{column_counter}=strcat('LC_MS fraction number_',Experimental_channel,'_',mat2str(column_counter));
end

fprintf(fid_summary_1, [Column_header, '\n'],...
'Experiment name', Column_name1{:}); %Write Header
fprintf(fid_summary_1, '%s,', Experiment_name1);
fprintf(fid_summary_1,'%6.4f,', re_calibrated_fraction_number(1,1:end));
fprintf(fid_summary_1,'\n');
fclose(fid_summary_1);

%Write out adjusted Raw data files for experimental channel tested and HvsM

%Experimental channel information
for column_counter=1:lenght_new_fraction
    Column_name{column_counter}=strcat('Ratio_',Experimental_channel,'_',mat2str(New_Fraction_numbering(column_counter)));
end

fid9_Name = strcat('Adjusted_',Experimental_channel,'_Raw_data_maxquant_rep',mat2str(alignment_counter),'.csv');
Column_header=repmat('%s,', 1, ((lenght_new_fraction)+2));

fid9 = fopen(fid9_Name,'at');
fprintf (fid9, [Column_header, '\n'],...
'Protein name', 'Replicate', Column_name{:});  %Write Header
for counter_fifteen=1:Raw_data_for_adjustment_size(1)
    fprintf(fid9, '%s,', Raw_data_import_for_adjustment.textdata{counter_fifteen+1,1});
    fprintf(fid9,'%6.4f,', alignment_counter); 
    fprintf(fid9,'%6.4f,', adjusted_raw_data(counter_fifteen,1:lenght_new_fraction));
    fprintf(fid9,'\n');
end
fclose(fid9);

%Write adjust table for Roc analysis

locations_of_protein_data=zeros(1,1);
%find location of protein to write out
for Roc_counter=1:length(Gaus_import_for_adjustment.textdata)-1
%find location of protein to write out
[location_Protein_in_Raw]=ind2sub(length(Raw_data_import_for_adjustment.textdata(:,1)),...
strmatch(Gaus_import_for_adjustment.textdata{Roc_counter+1,4},Raw_data_import_for_adjustment.textdata(:,1),'exact'));
locations_of_protein_data(Roc_counter)=location_Protein_in_Raw;
end


fid9B_Name = strcat('Adjusted_',Experimental_channel,'_Raw_for_ROC_analysis_rep',mat2str(alignment_counter),'.csv');
Column_header=repmat('%s,', 1, ((lenght_new_fraction)+2));

fid9B = fopen(fid9B_Name,'at');
fprintf (fid9B, [Column_header, '\n'],...
'Protein name', 'Replicate', Column_name{:});  %Write Header
for Roc_counter2 =1:length(locations_of_protein_data)
    fprintf(fid9B, '%s,', Raw_data_import_for_adjustment.textdata{locations_of_protein_data(Roc_counter2),1});
    fprintf(fid9B,'%6.4f,',alignment_counter); 
    fprintf(fid9B,'%6.4f,', adjusted_raw_data(locations_of_protein_data(Roc_counter2)-1,1:lenght_new_fraction));
    fprintf(fid9B,'\n');
end
fclose(fid9B);

%HvsM channel information Experimental channel
for column_counter=1:lenght_new_fraction
    Column_name{column_counter}=strcat('Ratio_HvsM_',mat2str(New_Fraction_numbering(column_counter)));
end

fid9B_Name = strcat('Adjusted_HvsM_Raw_data_maxquant_rep',mat2str(alignment_counter),'.csv');
Column_header=repmat('%s,', 1, ((lenght_new_fraction)+2));

fid9B = fopen(fid9B_Name,'at');
fprintf (fid9B, [Column_header, '\n'],...
'Protein name', 'Replicate', Column_name{:});  %Write Header
for counter_fifteenB=1:Raw_data_for_adjustment_size(1)
    fprintf(fid9B, '%s,', HvsM_Raw_data_import_for_adjustment.textdata{counter_fifteenB+1,1});
    fprintf(fid9B,'%6.4f,', HvsM_Raw_data_import_for_adjustment.data(counter_fifteenB,1)); 
    fprintf(fid9B,'%6.4f,', HvsM_adjusted_raw_data(counter_fifteenB,1:lenght_new_fraction));
    fprintf(fid9B,'\n');
end
fclose(fid9B);

%Write out adjusted Chromatograms
fid6_Name = strcat('Adjusted_Chromatograms_vobose_rep',mat2str(alignment_counter),'_','.csv');
fid6B = fopen(fid6_Name,'at');
for counter_fifteen=1:Protein_names_gaussian_for_adjustment_length
    fprintf(fid6B, '%s,', Gaus_import_for_adjustment.textdata{counter_fifteen+1,4});
    fprintf(fid6B,'%6.4f,', Combining_Adjusted_Chromatograms(counter_fifteen,:));
    fprintf(fid6B,'\n');
end
fclose(fid6B);

%Write out table with altered center
fid7_Name = strcat('Adjusted_Combined_OutputGaus_rep',mat2str(alignment_counter),'.csv');
fid7 = fopen(fid7_Name,'at');
fprintf (fid7,'%s,%s,%s,%s,%s,%s,%s\n',...
'Protein name', 'Height', 'Center','Width','SSE','adjrsquare', 'Complex Size');  %Write Header
for counter_sixteenth=1:Protein_names_gaussian_for_adjustment_length
    fprintf(fid7, '%s,', Gaus_import_for_adjustment.textdata{counter_sixteenth+1,4});
    fprintf(fid7,'%6.4f,', Gaus_import_for_adjustment.data(counter_sixteenth,:));
    fprintf(fid7,'\n');
end
fclose(fid7);
cd('../');
end


%Combine adjusted gaussian values for further analysis
%Import title information from first replicate
Adjusted_Gaus_import_Name_title= strcat('Realignment\Adjusted_Combined_OutputGaus_rep',mat2str(1),'.csv');
Adjusted_Gaus_import_title=importdata(Adjusted_Gaus_import_Name_title);
Title_Gaus_import=regexp(Adjusted_Gaus_import_title.textdata{1,1}, ',' ,'split');

fid10_Name = strcat('Adjusted_',Experimental_channel,'_Combined_OutputGaus.csv');
fid10 = fopen(fid10_Name,'at');
fprintf (fid10,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n',...
'Guassian index number', 'Protein_number','Replicate',Title_Gaus_import{:});  %Write Header

for alignment_counter=1:Unique_replicate
    
%Import Matlab processed data- alignment data
Adjusted_Gaus_import_Name= strcat('Realignment\Adjusted_Combined_OutputGaus_rep',mat2str(alignment_counter),'.csv');
Adjusted_Gaus_import=importdata(Adjusted_Gaus_import_Name);           %Import the Gaus data from previous scripts
Adjusted_protein_names_gaussian_length=length(Adjusted_Gaus_import.data);

%Import Matlab processed data- header data
Gaus_import_Name= strcat('Processed Gaussian\',Experimental_channel,'_Combined_OutputGaus_rep',mat2str(alignment_counter),'.csv');
Gaus_import=importdata(Gaus_import_Name);           %Import the Gaus data from previous scripts

%Write to file
for combined_gaus_writeout=1:Adjusted_protein_names_gaussian_length
fprintf(fid10, '%6f,', str2num(Gaus_import.textdata{combined_gaus_writeout+1,1}));
fprintf(fid10, '%6f,', str2num(Gaus_import.textdata{combined_gaus_writeout+1,2}));
fprintf(fid10, '%6f,', str2num(Gaus_import.textdata{combined_gaus_writeout+1,3})); 
fprintf(fid10, '%s,', Adjusted_Gaus_import.textdata{combined_gaus_writeout+1,1});
fprintf(fid10,'%6.4f,', Adjusted_Gaus_import.data(combined_gaus_writeout,:));
fprintf(fid10,'\n');

end
end
fclose(fid10);

%Move files to folder for further analysis
movefile(fid10_Name,'Combined Realigned data');

%Combined adjusted raw data
%Import title information from first replicate
Adjusted_Raw_data_Name_title= strcat('Realignment\Adjusted_',Experimental_channel,'_Raw_data_maxquant_rep',mat2str(1),'.csv');
Adjusted_Raw_data_title=importdata(Adjusted_Raw_data_Name_title);

%create header information
Title_import2=Adjusted_Raw_data_title.textdata(1,:);
Column_header=repmat('%s,', 1, (length(Title_import2)));

fid11_Name = strcat('Adjusted_',Experimental_channel,'_Raw_data_maxquant.csv');
fid11 = fopen(fid11_Name,'at');
fprintf (fid11,[Column_header, '\n'], Title_import2{:});  %Write Header

for alignment_counter=1:Unique_replicate
    
%Import Matlab processed data- alignment data
Adjusted_Raw_data_import_Name= strcat('Realignment\Adjusted_',Experimental_channel,'_Raw_data_maxquant_rep',mat2str(alignment_counter),'.csv');

%import data files summary data
f5 = fopen (Adjusted_Raw_data_import_Name);
Adjusted_Raw_data = textscan(f5, '%s', 'Delimiter',',');
fclose(f5);

%Reshape data for use in analysis
Adjusted_Raw_Dimension=size(Adjusted_Raw_data{:});
Adjusted_Raw_data_reshaped=reshape(Adjusted_Raw_data{:},68,(Adjusted_Raw_Dimension(1)/68))';
Adjusted_Raw_data_reshaped_size=size(Adjusted_Raw_data_reshaped);

%Write to file
for combined_gaus_writeout=2:Adjusted_Raw_data_reshaped_size(1)
fprintf(fid11, '%s,',  Adjusted_Raw_data_reshaped{combined_gaus_writeout,1});
fprintf(fid11, '%6f,', str2num(Adjusted_Raw_data_reshaped{combined_gaus_writeout,2}));
fprintf(fid11, '%s,',  Adjusted_Raw_data_reshaped{combined_gaus_writeout,3:end}); 
fprintf(fid11,'\n');

end
end
fclose(fid11);

%Move files to folder for further analysis
movefile(fid11_Name,'Combined Realigned data');


%Combined adjusted raw data
%Import title information from first replicate
Adjusted_Raw_data_Name_title= strcat('Realignment\Adjusted_HvsM_Raw_data_maxquant_rep',mat2str(1),'.csv');
Adjusted_Raw_data_title=importdata(Adjusted_Raw_data_Name_title);

%create header information
Title_import2=Adjusted_Raw_data_title.textdata(1,:);
Column_header=repmat('%s,', 1, (length(Title_import2)));

fid11B_Name = strcat('Adjusted_HvsM_Raw_data_maxquant.csv');
fid11B = fopen(fid11B_Name,'at');
fprintf (fid11B,[Column_header, '\n'], Title_import2{:});  %Write Header

for alignment_counter=1:Unique_replicate
    
%Import Matlab processed data- alignment data
Adjusted_Raw_data_import_Name= strcat('Realignment\Adjusted_HvsM_Raw_data_maxquant_rep',mat2str(alignment_counter),'.csv');

%import data files summary data
f5 = fopen (Adjusted_Raw_data_import_Name);
Adjusted_Raw_data = textscan(f5, '%s', 'Delimiter',',');
fclose(f5);

%Reshape data for use in analysis
Adjusted_Raw_Dimension=size(Adjusted_Raw_data{:});
Adjusted_Raw_data_reshaped=reshape(Adjusted_Raw_data{:},68,(Adjusted_Raw_Dimension(1)/68))';
Adjusted_Raw_data_reshaped_size=size(Adjusted_Raw_data_reshaped);

%Write to file
for combined_gaus_writeout=2:Adjusted_Raw_data_reshaped_size(1)
fprintf(fid11B, '%s,',  Adjusted_Raw_data_reshaped{combined_gaus_writeout,1});
fprintf(fid11B, '%6f,', str2num(Adjusted_Raw_data_reshaped{combined_gaus_writeout,2}));
fprintf(fid11B, '%s,',  Adjusted_Raw_data_reshaped{combined_gaus_writeout,3:end}); 
fprintf(fid11B,'\n');

end
end
fclose(fid11B);

%Move files to folder for further analysis
movefile(fid11B_Name,'Combined Realigned data');

%Combined adjusted fraction data
%Import title information from first replicate
Adjusted_fraction_Name_title= strcat('Realignment\fraction_overlap summary_',Experimental_channel,'_replicate_',mat2str(1),'.csv');
Adjusted_fraction_title=importdata(Adjusted_fraction_Name_title);

%create header information
Title_import3=Adjusted_fraction_title.textdata(1,:);
Column_header=repmat('%s,', 1, (length(Title_import3)));

fid12_Name = strcat('Adjusted_fraction_overlap summary_',Experimental_channel,'.csv');
fid12 = fopen(fid12_Name,'at');
fprintf (fid12,[Column_header, '\n'], Title_import3{:});  %Write Header

for alignment_counter=1:Unique_replicate
    
%Import Matlab processed data- alignment data
Adjusted_fraction_import_Name= strcat('Realignment\fraction_overlap summary_',Experimental_channel,'_replicate_',mat2str(alignment_counter),'.csv');

%import data files summary data
f6 = fopen (Adjusted_fraction_import_Name);
Adjusted_fraction_import = textscan(f6, '%s', 'Delimiter',',');
fclose(f6);

%Reshape data for use in analysis
Adjusted_fraction_import_Dimension=size(Adjusted_fraction_import{:});
%Determine the dimension of the array
reshap_values=Adjusted_fraction_import_Dimension(1)/2 %divide by two as have two rows
Adjusted_fraction_import_reshaped=reshape(Adjusted_fraction_import{:},reshap_values,(Adjusted_fraction_import_Dimension(1)/reshap_values))';
Adjusted_fraction_import_reshaped_size=size(Adjusted_fraction_import_reshaped);

%Write to file
for combined_gaus_writeout=2:Adjusted_fraction_import_reshaped_size(1)
fprintf(fid12, '%s,',  Adjusted_fraction_import_reshaped{combined_gaus_writeout,1});
for combined_gaus_writeout2=2:Adjusted_fraction_import_reshaped_size(2)
fprintf(fid12, '%6f,', str2num(Adjusted_fraction_import_reshaped{combined_gaus_writeout,combined_gaus_writeout2}));
end
fprintf(fid12,'\n');

end
end
fclose(fid12);

%Move files to folder for further analysis
movefile(fid12_Name,'Combined Realigned data');
fclose(fid6C);
fclose(fid6D);

cd('../')
end


fclose('all');