%%%%% 17/10/11
%%%% Thermal Shift Evaluation Software - MTSA
%%% Michele N. Schulz
%%% michele@ysbl.york.ac.uk

%%% Please cite:
%%% "MTSA - A Matlab Programme to Fit Thermal Shift Data"
%%% Michele N. Schulz, Jens Landstrom, Roderick E. Hubbard
%%% Analytical Biochemistry, 2011

% Copyright 2011 Michele Schulz
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 
% Unless required by applicable law or agreed to in writing, software distributed 
% under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR 
% CONDITIONS OF ANY KIND, either express or implied. See the License for the 
% specific language governing permissions and limitations under the License. 




function varargout = MTSA_111111_Boltzmann(varargin)
% MTSA_111111_Boltzmann M-file for MTSA_111111_Boltzmann.fig
%      MTSA_111111_Boltzmann, by itself, creates a new MTSA_111111_Boltzmann or raises the existing
%      singleton*.
%
%      H = MTSA_111111_Boltzmann returns the handle to a new MTSA_111111_Boltzmann or the handle to
%      the existing singleton*.
%
%      MTSA_111111_Boltzmann('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MTSA_111111_Boltzmann.M with the given input arguments.
%
%      MTSA_111111_Boltzmann('Property','Value',...) creates a new MTSA_111111_Boltzmann or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before MTSA_111111_Boltzmann_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to MTSA_111111_Boltzmann_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help MTSA_111111_Boltzmann

% Last Modified by GUIDE v2.5 29-Nov-2011 16:13:35

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @MTSA_111111_Boltzmann_OpeningFcn, ...
                   'gui_OutputFcn',  @MTSA_111111_Boltzmann_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before MTSA_111111_Boltzmann is made visible.
function MTSA_111111_Boltzmann_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to MTSA_111111_Boltzmann (see VARARGIN)

% Choose default command line output for MTSA_111111_Boltzmann
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes MTSA_111111_Boltzmann wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = MTSA_111111_Boltzmann_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;



%%%%%Table to enter the fluorescence data


% --- Executes when selected cell(s) is changed in uitable1.
function uitable1_CellSelectionCallback(hObject, eventdata, handles)
% hObject    handle to uitable1 (see GCBO)
% eventdata  structure with the following fields (see UITABLE)
%	Indices: row and column indices of the cell(s) currently selecteds
% handles    structure with handles and user data (see GUIDATA)


% --- Executes when entered data in editable cell(s) in uitable1.
function uitable1_CellEditCallback(hObject, eventdata, handles)
% hObject    handle to uitable1 (see GCBO)
% eventdata  structure with the following fields (see UITABLE)
%	Indices: row and column indices of the cell(s) edited
%	PreviousData: previous data for the cell(s) edited
%	EditData: string(s) entered by the user
%	NewData: EditData or its converted form set on the Data property. Empty if Data was not changed
%	Error: error string when failed to convert EditData to appropriate value for Data
% handles    structure with handles and user data (see GUIDATA)



%%%%Paste Data menu 

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

import =importdata('-pastespecial');
set(handles.uitable1,'Data',import);



%%%%%%%%%%%How many wells to fit?

function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double


% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end





%%%%%%%%%%%%%%%Which well number is first blank:

function edit7_Callback(hObject, eventdata, handles)
% hObject    handle to edit7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit7 as text
%        str2double(get(hObject,'String')) returns contents of edit7 as a double


% --- Executes during object creation, after setting all properties.
function edit7_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end





%%%%%%%%%%%%%%Spacing between blanks

function edit3_Callback(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double


% --- Executes during object creation, after setting all properties.
function edit3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


%%%%%%%%%%%Starting temperature

function edit4_Callback(hObject, eventdata, handles)
% hObject    handle to edit4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit4 as text
%        str2double(get(hObject,'String')) returns contents of edit4 as a double


% --- Executes during object creation, after setting all properties.
function edit4_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


%%%%%%%%Increment of temperature

function edit5_Callback(hObject, eventdata, handles)
% hObject    handle to edit5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit5 as text
%        str2double(get(hObject,'String')) returns contents of edit5 as a double


% --- Executes during object creation, after setting all properties.
function edit5_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



%%%%%%%Number of cycles on QPCR while increasing the temperature

function edit6_Callback(hObject, eventdata, handles)
% hObject    handle to edit6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit6 as text
%        str2double(get(hObject,'String')) returns contents of edit6 as a double


% --- Executes during object creation, after setting all properties.
function edit6_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




function edit9_Callback(hObject, eventdata, handles)
% hObject    handle to edit9 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit9 as text
%        str2double(get(hObject,'String')) returns contents of edit9 as a double


% --- Executes during object creation, after setting all properties.
function edit9_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit9 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end





%%%%Click to fit button

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)



%%%Script starts here



%%%clear command window

clc


%%%%%%%%%%%%%%%%GET USER ENTRIES

%%Get the input table
inputTable = get(handles.uitable1, 'Data');



%%How many first curves to do (max for a plate is 96 -  obviously)
n = str2double(get(handles.edit1,'String'));



%%Repeated Blank
blankRepeat = str2double(get(handles.edit3,'String'));

%%First occurence of blank
firstBlank = str2double(get(handles.edit7,'String'));


%%get temperature data
startTemp = str2double(get(handles.edit4,'String')); 
increment = str2double(get(handles.edit5,'String')); 
cycles = str2double(get(handles.edit6,'String')); 


%%Filename
MyResults = get(handles.edit9,'String');


%%%%%%%%%%%%%%%%%USER ENTRIES END

%%delete old diary file if present to not append to old data
if exist(MyResults) == 2
    delete(MyResults);
end



%%%%write all output into a fresh diary matlab file
diary(MyResults)   



%%%%%%%%%%%
%all is a matrix of data from input table
all = [inputTable];

%%create temperature vector
endTemp = startTemp + cycles -1;

%temperature vector, transposed
xdata = [startTemp:increment:endTemp]';


%%expected melting temperature
Texpected = 60; 



%%clear old plot if stil open
clf


%%to make the output nicer print that errors are following
disp('Occuring errors below:');
disp(' ')

%how to plot the original data. in this case blue dots
DataLineSpec = '.b';

%indicates how the fit should be shown, here magenta dots
FitLineSpec = '.m';


%%start loop over all data

for i = 1:n
    
 
    %%%clear plot after each fit

    cla %%clear contents of plot

  %%Count variable
Counter = num2str(i);    
    
%%%try to catch the error    
try

    %%extract data
    ydata = all(:,i);
    

    %Boltzmann fit
    Boltzmann = fittype('min + (max-min)./(1+exp((Tm-x)/a))');
    

    
    %%define the area of the curve which should be fitted:
    
    %%find first the maximum
    maxydata = max(ydata);
    %%find index of maximum and define it as border
    indexRightBorder = find(ydata == maxydata);
    
    %%%just in case there is more than one element with the same value... just
    %%%take the first one        
    indexRightBorder = indexRightBorder(1,:);
    
    
    %%%Look only at the area left of maximum
    Leftydata = ydata([1:indexRightBorder]);
    
    %%%find the minimum left of the maximum
    minydata = min(Leftydata);    
    
    %%%find the index of the minimum
    indexLeftBorder = find(minydata == ydata);    
    %%%just in case there is more than one element with the same value... just
    %%%take the first one
    indexLeftBorder = indexLeftBorder(1,:);
    
    
    
    %%%get the corresponding values on the x-axis
    leftBorder = xdata(indexLeftBorder);
    rightBorder = xdata(indexRightBorder);


    
    %%use these values to define the borders of the to excluded data
    excludedData = excludedata(xdata,ydata,'domain',[leftBorder rightBorder]);
    
    
    %create fit options structure
    myoptions = fitoptions('Exclude', excludedData,...
        'Method','NonlinearLeastSquares',...
        'Lower',[Texpected - 50,0,0.98*maxydata,0.98*minydata],...
        'Upper',[Texpected + 50,10,1.02*maxydata,1.02*minydata],...
        'StartPoint',[Texpected,5,maxydata,minydata]);
    
    
    
    %%Fit and calculate error statistics and gives output
    [FittedCurve,Goodness,Output] = fit(xdata,ydata,Boltzmann,myoptions);
    
    %%assign variables which change with the loop
    eval([ 'FittedCurve' num2str(Counter) ' = FittedCurve;' ]);
    eval([ 'Goodness' num2str(Counter) ' = Goodness;' ]);
    eval([ 'Output' num2str(Counter) ' = Output;' ]);
    
    
    %plot fit and datapoints in defined colours and label
    hold on %%%plots in same window as before
    p = plot(FittedCurve, FitLineSpec, xdata, ydata);
    
    xlabel('Temperature','FontSize',14);
    ylabel('Fluorescence','FontSize',14);
    title(['Thermal Melting Curve of Well ' num2str(Counter)],'FontSize',14)
    h_legend = legend('Experimental Data','Fitted Data','Location','NorthWest');
    set(h_legend,'FontSize',10);
   
    
    %%create string for plot

    string1 = 'Tm=';
    string2 = eval(['FittedCurve' num2str(Counter) '.Tm']);  %%% Tm =T m
    string3 = num2str(string2);
    
    string4 = 'Rsquare=';
    string5 = eval(['Goodness' num2str(Counter) '.rsquare']);
    string6 = num2str(string5);
    
    stringTm = strcat(string1,string3);
    stringRsquare = strcat(string4,string6);
    
    %%print the text on the plot
    text(endTemp-10, 0.95*maxydata,{stringTm;stringRsquare},...  
       'VerticalAlignment','middle',...
       'HorizontalAlignment','left',...        
       'FontSize',10);
    
    
        
    
    
    %%%Save as *.png
    %%Give curve a name with number, e.g. 1 to 96 for one plate
    FileNamePng = strcat('FittedCurve','_',Counter,'.png');
    
    X=getframe(gcf);
    isempty(X.colormap);
    imwrite(X.cdata, FileNamePng);
    
    
    
%%%catch the error    
catch exception
disp(['Curve number ' Counter ' could not be fitted. Check data!']); 
end


end





%%To make output nicer, indicate what's coming next
disp('All fits with statistics:');
disp(' ')



%%output only the relevant values
%%catch error when it does not work
for j = 1:n
        try  exist(eval (['FittedCurve' num2str(j)]))   == 1;
        eval( ['FittedCurve' num2str(j)]);
        eval(['Goodness' num2str(j)]);
        eval(['Output' num2str(j)]);
        catch exception
        disp('No fit available'); 
    end
end


%%output only the Tms
%%catch error when it does not work
for j = 1:n
        try  exist(eval (['FittedCurve' num2str(j)]))   == 1;
    
      struct(['Tm' num2str(j)] ,...    
        eval(['FittedCurve' num2str(j) '.Tm']),...  %%% Tm =Tm
        ['Rsquare' num2str(j)] ,eval( ['Goodness' num2str(j) '.rsquare']) )
       
        catch exception
        disp('No fit available'); 
    end
end


%%%show only the Tms 

for j = 1:n
        try  exist(eval (['FittedCurve' num2str(j)]))   == 1;   
     disp( struct(['Tm' num2str(j)] ,  eval(['FittedCurve' num2str(j) '.Tm']) ))  %%% Tm =Tm
       catch exception
        disp('No fit available'); 
    end
end




%%%Get the Average Tm for the blanks

 TmSum = 0;
     
%%catch the error for the case that some blanks could not be fitted   
     
        for i=1:n;
          try
 
             if mod(i,blankRepeat) == firstBlank
                 TmSum = TmSum + eval(['FittedCurve' num2str(i) '.Tm']);  %%% Tm =Tm
             end
         
          catch exception
             disp('Check your blanks');
          end
         
        end
         
%%the avergae of all blanks     
AverageTmOfBlanks = TmSum / (n / blankRepeat)






%%%And the standard deviation of the blanks is:
 BlankTms = [];
 
     for i=1:n;
         if mod(i,blankRepeat) == firstBlank
            BlankTms = [BlankTms, eval(['FittedCurve' num2str(i) '.Tm'])];  %%% Tm =Tm
         end 
     end
     
     
 StandardDeviationOfBlank = std(BlankTms)

 
 %%Give a warning when the standard deviation of the blanks is over 0.3
 
 if StandardDeviationOfBlank > 0.3
     disp('WARNING: The standard deviation is higher than 0.3!');
 end
 
 disp(' ')
 
 
 %%%%show Tms with thermal shift in respect to the blanks and the Rsquare
 %warning if values are not good

disp('Final fitting summary below in the format');
disp('FitNumber: [Tm, shift in respect to blank, Rsquare:');
disp(' ')


 for j = 1:n
     try  exist(eval (['FittedCurve' num2str(j)]))   == 1;   
         
     disp(struct(  ['Fit' num2str(j)] , [eval(['FittedCurve' num2str(j) '.Tm']),... %%% Tm =Tm
        eval(['FittedCurve' num2str(j) '.Tm'])- AverageTmOfBlanks,... %%% Tm =Tm
        eval( ['Goodness' num2str(j) '.rsquare'] ) ]   )  ); 
  
     if eval( ['Goodness' num2str(j) '.rsquare'] ) < 0.999
        disp(['WARNING: Fit' num2str(j) ' IS NOT WELL FITTED'])
        disp(' ')    
     end  
    
     catch exception
     disp( struct(  ['Tm' num2str(j)] , 'No fit available' ) ); 
     end
 end 
 

 
%%%%close diary file 
diary off




%%delete the redundant loop variables
clear ('FittedCurve', 'Goodness','Output','xdata',',ydata,',...
    'maxydata','minydata','X','Title','excludedData','myoptions',...
    'Counter','DataLineSpec','FileNamePng','FitLineSpec',...
    'indexLeftBorder','indexRightBorder','p','leftBorder',...
    'rightBorder','ydata','i''n','j')
