# -*- coding: utf-8 -*- """hrm_analysis.ipynb Automatically generated by Colaboratory. Original file is located at https://colab.research.google.com/drive/1vF94CWORVtXSZDjr0zFpaoW-O2YvsZyO """ import pandas as pd import numpy as np import matplotlib.pyplot as plt from scipy.signal import savgol_filter import math from numpy.polynomial.polynomial import Polynomial data = pd.read_csv("HRM.csv", index_col=0, header=0, sep= ";", decimal=",") #use the separator according to your file. Common separator are ",", ";" data.loc[:,'degree']=data.index data.index = pd.to_datetime(data.index, unit='ms') folder_results = "results" # change this to your results folder to save figures def plot_curve(index, N_sample_1, B_T, sample_1, name, peak): plt.plot(index, N_sample_1, label='normalized_raw') plt.plot(index, N_M, label='normalized_EBS') plt.plot(index, -1*B_T, label='Background') plt.plot(index, M_T, label='Raw MT') plt.plot(index, sample_1, label='Raw sample') plt.legend() plt.suptitle(name) plt.title("Melting temp" + repr(index[peak[-1]])) #Peaks on the derivative curves can be used for the identification of the melting temperatures (Tm) plt.savefig(results_folder+"/sample_"+name+"_melt_curve") plt.show() def plot_derivative(index, derivative_values, key, peak): ax = plt.plot(index, derivative_values, c='red') plt.title(key + ": Negative First Derivative Curve") plt.axvline(start, linestyle='--', c='blue') plt.axvline(end, linestyle='--', c='black') for p in peak: plt.axhline(derivative_values[p], linestyle='--', c='grey') plt.text(index[p], derivative_values[p], "temp: "+"{:.2f}".format(index[p])) plt.ylabel("-dF/dT") plt.xlabel("Temperature") plt.savefig(folder+"/"+key+"_derivative") plt.show() def plot_melting_region(index, sample_1, key, start, end): plt.plot(index, sample_1, c='red') plt.title(key) plt.axvline(start, linestyle='--', c='blue') plt.axvline(end, linestyle='--', c='black') plt.xlabel("Temperature") plt.ylabel("Fluorescence") plt.text(start, np.mean(sample_1)+0.2, "T_m_start $"+"{:.2f}".format(start)+"$") plt.text(end, np.mean(sample_1)-0.1, "T_m_end $"+"{:.2f}".format(end)+"$") plt.savefig(folder_results+"/"+key) plt.show() """**Defining samples and negative control** Following the example, define which columns belong to each sample and which columns are the negative control """ samples = dict() samples['Control_1'] = ['A2 ', 'A3 '] #Sample "Control_1" has two replicas in columns A2 and A3 samples['Control_2'] = ['A4 ', 'A5 '] #Sample "Control_2" has two replicas in columns A4 and A5 samples['A'] = ['A6 ', 'A7', 'A8'] #Sample "A" has three replicas in columns A6, A7 and A8 samples['Negative'] = ['C2 ','C3 '] #Negative control is in columns C2 and C3 control_samples= ['Control_1', "Control_2"] #Add your control samples name as you named them before. Separate them using , """**HRM Analysis** In this part, we process each sample using the algorithm described in the paper. The three steps are signaled by comments: finding melting region, apply EBS and Melting peak identification """ all_curves = pd.DataFrame() all_derivatives = pd.DataFrame() for key in samples: try: columns = list(filter(lambda x: any(s in x for s in samples[key]), data.columns)) sample = data[columns] sample.loc[:,'degree']=data['degree'] sample = sample.resample('1us').mean().interpolate() sample = sample.resample('60us').ffill() sample = sample.dropna() #find RecWindow index = sample['degree'].values rec_window = len(index[index<index[0]+1]) #number of samples within 1 degree diff = index - np.roll(index, 1) delta = np.min(diff[1:]) sample_1 = np.median(sample, axis=1) #we use median of replicas instead of mean derivative_values = -1*savgol_filter(sample_1, delta = delta, window_length=rec_window*2+1, polyorder=2,deriv=1) #polynomial order is set to 2, window is number of samples within 1 degree RecWindow*2+1 #Finding melting region tangent_derivative = savgol_filter(derivative_values,delta= delta/100, window_length=rec_window*2+1, polyorder=2,deriv=1) #polynomial order is set to 2, window is number of samples within 1 degree RecWindow*2+1 angles = np.degrees(np.arctan(tangent_derivative)) start = index[(angles>50)&(angles<55)] end = index[(angles<-30)&(angles<35)] if len(start)>0 and len(end)>0: #at least one melting region was found start = start[0] end = end[end>start] end = end[-1] #plot melting region plot_melting_region(index, sample_1, key, start, end) #apply EBS for melting curve #the raw fluorescence F(T) is modeled as composing of the melting curve M(T) and an exponentially decaying background B(T). It can be expressed using the following mathematical model #F(T) = M(T) + B(T) #first determine TR and Tl where derivative of M(T) is 0, outside of start and end TL = start - 2 TR = end + 2 tl_i = 0 tr_i = len(index)-1 if len(np.where(index<TL)[0])>0: tl_i = np.where(index<TL)[-1][-1] else: TL = index[0] if len(np.where(index>TR)[0])>0: tr_i = np.where(index>TR)[0][0] else: TR = index[-1] #We remove data that is not in the melting region sample_1 = sample_1[tl_i:tr_i+1] n_index = index[tl_i:tr_i+1] #Melting peak identification peak = np.where(np.sign(tangent_derivative[:-1]) > np.sign(tangent_derivative[1:]))[0] #tm is where there is a change from positive to negative in tangent derivative tangent_derivative = tangent_derivative[tl_i:tr_i+1] plot_derivative(index, derivative_values, key, peak) index = n_index #Get Background a = (math.log(derivative_values[tr_i])-math.log(derivative_values[tl_i]))/(TR-TL) C = derivative_values[tl_i] / a exponents = a*(index-TL) B_T =C*np.exp(exponents) M_T = sample_1 - B_T N_M = 100 * (M_T-np.min(M_T))/(np.max(M_T)-np.min(M_T)) #Normalize sample N_sample_1 = 100 * (sample_1-np.min(sample_1))/(np.max(sample_1)-np.min(sample_1)) #Plot samples and normalized plot_curve(index, N_sample_1, B_T, sample_1, key, peak) #save the curve for plotting all together sample_data = pd.DataFrame(N_M, index=index, columns=[key]) all_curves = pd.concat([all_curves, sample_data], axis=1) derivative = -1*savgol_filter(N_M, delta=delta, window_length=rec_window*2+1, polyorder=2,deriv=1) sample_data = pd.DataFrame(derivative, index=index, columns=[key]) all_derivatives = pd.concat([all_derivatives, sample_data], axis=1) except (Exception) as e: print("error on sample " + key) print(e) sample_data = pd.DataFrame(np.ones(len(index)), index=index, columns=[key]) all_derivatives = pd.concat([all_derivatives, sample_data], axis=1) continue """**Plot results all together**""" #plot difference plot import seaborn as sns reshaped_df = all_curves.reset_index() reshaped_df = pd.melt(reshaped_df, id_vars='index') g= sns.lineplot(x='index', y='value', hue='variable', data=reshaped_df) plt.title("All samples melt curve") plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) g.get_figure().savefig(results_folder+"/melting_all.png", bbox_inches='tight', dpi=300) plt.show() reshaped_df = all_derivatives.reset_index() reshaped_df = pd.melt(reshaped_df, id_vars='index') g= sns.lineplot(x='index', y='value', hue='variable', data=reshaped_df) plt.title("All samples negative derviative curve") plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) g.get_figure().savefig(results_folder+"/neg_derivative_all.png",, bbox_inches='tight', dpi=300) plt.show() med_curve = np.mean(all_curves[control_samples], axis=1) plt.plot(med_curve) plt.show() all_curves_sub =all_curves.sub(med_curve, axis='index') reshaped_df = all_curves_sub.reset_index() reshaped_df = pd.melt(reshaped_df, id_vars='index') g= sns.lineplot(x='index', y='value', hue='variable', data=reshaped_df) plt.title("Difference plot") plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) g.get_figure().savefig(results_folder+"/difference_melt.png", bbox_inches='tight', dpi=300) plt.show() #Difference plot with all temperatures med_curve = med_curve.fillna(0) all_curves_sub =all_curves.sub(med_curve, axis='index') all_curves_sub = all_curves_sub.iloc[20:,] reshaped_df = all_curves_sub.reset_index() reshaped_df = pd.melt(reshaped_df, id_vars='index') g= sns.lineplot(x='index', y='value', hue='variable', data=reshaped_df) plt.title("Difference plot with all temperatures") plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) g.get_figure().savefig(results_folder+'/difference_melt_all_temps.png', bbox_inches='tight', dpi=300) plt.show()