%%Analyze Calcium Influx Data %% By Yingjie Du %%Identify Peaks and valleys from Graphs %%Isolate those data points %%Export those peaks and values to get average max intensity and amplitude %%Obtain frequency by calculating number of peaks %% Variables that need to be changed: %% z: creates sections to identify lowest baseline (should be relatively small) %% z: creates sections for identifying highest data points %% cutoff multiplier: Identifies data points/peaks above a certain threshold %% 3 Import Titles, 1 Export title and 1 Figure title %% Use this section for isolating peak/valley points %Load data and separate sets of data load Example_mSMCs_WT_04.txt time_raw = Example_mSMCs_WT_04(:,1); %Everything in first column and all rows max_intensity_raw = Example_mSMCsWT04(:,2); %Everything in second column and all rows time = time_raw'; intensity = max_intensity_raw'; plot(time,intensity); hold on s = size(intensity); s = s(1,2); %%Filter 1: Eliminate background noise to find peaks that matter %%2 Things that can be changed: increment, cutoff multiplier Condition: z %%has to be factor of new_max_intensity size %% This section will analyze entire curve for peaks %% Calculations: Frequency g = 1; h = 1; m = 1; x = 1; z = 2; %increment for baseline (choose healthy sample size!!!!) for L = 1:z:s-1 local_average_intensity(1,g) = mean(intensity(1,L:L+z-1)); g = g+1; end %global_average_intensity = min(local_average_intensity); global_average_intensity = 556; baseline_estimate = .03*global_average_intensity;%intensity cutoff for intensity change calculation only %want smaller baseline estimate z = 50; % increment for peaks (tune as needed)!!!!!! for L = 1:z:s-1 for i = 1:z if i+L-1 > 1 && i+L-1 < s %peak not the last or first data point if intensity(1,i+L-1) > 1.255*global_average_intensity %intensity cutoff for peak identification only if intensity(1,i+L-1) > intensity(1,i+L) && intensity(1,i+L-1) > intensity(1,i+L-2) %isolate peak final_intensity(1,m) = intensity(1,i+L-1); %Save highest peaks final_time(1,m) = time(1,i+L-1); %save points left and right of peak for future %calculations peak_intensity2(1,h+1) = intensity(1,i+L-1); peak_time2(1,h+1) = time(1,i+L-1); o = i; q = i; if o+L-2 == 0 o = o+1; end check_diff1 = abs(intensity(1,o+L-2)-global_average_intensity); while check_diff1 > baseline_estimate && o+L-2 < s o = o-1; check_diff1 = abs(intensity(1,o+L-2)-global_average_intensity); end peak_intensity2(1,h) = intensity(1,o+L-2); peak_time2(1,h) = time(1,o+L-2); check_diff2 = abs(intensity(1,q+L)-global_average_intensity); while check_diff2 > baseline_estimate && q+L < s q = q+1; check_diff2 = abs(intensity(1,q+L)-global_average_intensity); end peak_intensity2(1,h+2) = intensity(1,q+L); peak_time2(1,h+2) = time(1,q+L); h = h+3; m = m+1; end end end end end plot(final_time,final_intensity) hold on plot(peak_time2,peak_intensity2) hold on p = size(time_raw'); p = p(1,2); %array for global average value for k = 1:p A(1,k) = global_average_intensity; end plot(time_raw,A) hold on %% This next section will analyze each peak in detail based with previously %% made arrays %% Calculations: average intensity changes,time between peaks,peak width %% time b=1; for n = 2:h-1 if peak_intensity2(1,n) > peak_intensity2(1,n-1) && peak_intensity2(1,n) > peak_intensity2(1,n+1) %Calculation: intensity change peak_base = .5*(peak_intensity2(1,n-1) + peak_intensity2(1,n+1)); intensity_change(1,b) = peak_intensity2(1,n) - peak_base; %Calculation: time of peak width time_of_peak(1,b) = peak_time2(1,n+1) - peak_time2(1,n-1); %Calculation: average intensity change if n < h-3 time_between_peak(1,b) = peak_time2(1,n+2) - peak_time2(1,n+1); end b = b+1; end end if b == 2 time_between_peak (1,1) = 0; end %Calculating average intensity and stdev average_max_intensity = mean(final_intensity); stdev_max_intensity = std(final_intensity); %% Use this section for any remaining calculations and exporting data title('mSMCs_WT_04, cell 3') xlabel('Time(ms)') ylabel('Intensity') %calculating frequency = defined as the number of "flashes" we get within %time period. final_s = size(intensity_change); final_s = final_s(1,2); frequency = (final_s/60); %Calculating area under the curve andaverage intensity change based on %single baseline r_size = size(peak_intensity2); r = r_size(1,2); v = 1; c = 1; for d = 2:r-1 if peak_intensity2(1,d) > peak_intensity2(1,d-1)&& peak_intensity2(1,d)> peak_intensity2(1,d+1) base = peak_time2(1,d+1) - peak_time2(1,d-1); height = peak_intensity2(1,d) - .5*(peak_intensity2(1,d+1) + peak_intensity2(1,d-1)); area_curve(1,v) = .5*base*height; v = v+1; alt_int_change(1,c) = peak_intensity2(1,d)- global_average_intensity; c = c+1; end end % export_data = [average_max_intensity, stdev_max_intensity, max_intensity, min_intensity, difference_final_intensity, frequency]; export_data1 = final_intensity'; export_data2 = intensity_change'; export_data3 = alt_int_change'; export_data4 = time_of_peak'; export_data5 = time_between_peak'; export_data6 = area_curve'; export_data7 = frequency'; export_data8 = peak_intensity2'; export_data9 = peak_time2'; export_data10 = time'; export_data11 = intensity'; xlswrite('export mSMCs_WT_04.xls', export_data1, 'Sheet 1', 'A1') xlswrite('export mSMCs_WT_04.xls', export_data2, 'Sheet 1', 'B1') xlswrite('export mSMCs_WT_04.xls', export_data3, 'Sheet 1', 'C1') xlswrite('export mSMCs_WT_04.xls', export_data4, 'Sheet 1', 'D1') xlswrite('export mSMCs_WT_04.xls', export_data5, 'Sheet 1', 'E1') xlswrite('export mSMCs_WT_04.xls', export_data6, 'Sheet 1', 'F1') xlswrite('export mSMCs_WT_04.xls', export_data7, 'Sheet 1', 'G1') xlswrite('export mSMCs_WT_04.xls', export_data8, 'Sheet 1', 'H1') xlswrite('export mSMCs_WT_04.xls', export_data9, 'Sheet 1', 'I1') xlswrite('export mSMCs_WT_04.xls', export_data10, 'Sheet 1', 'J1') xlswrite('export mSMCs_WT_04.xls', export_data11, 'Sheet 1', 'K1')