%% To do %------------------------------------------------------------------- % enter warnings for: % attempting to analyze input-output .abf files % attempting to analyze noisy data % attempting to analyze data that does not follow electrophysio assay % test fit curve(done 10/6/10) % create output of stArray into txt file(done 10/6/10) % create progress bar(done 10/6/10) % create plots of important data(done 10/7/10) %------------------------------------------------------------------- %% Initialize variables %------------------------------------------------------------------- %create a cell array for each sweep (150 sweeps collected in 10908001) %and only plot the first 7 ms then find markers clear all; close all; clc; %initialize progress bar h = waitbar(0,'Initializing Global waitbar...',... 'position',[400 100 275 50]); button = []; button1 = []; button2 = []; button3 = []; button4 = []; button5 = []; button6 = []; %create a directory of all abf files directory = dir('*abf'); %initialize counters for figure plots and waitbar fig = 0; step = 0; %------------------------------------------------------------------- %% Load data %------------------------------------------------------------------- %read in each experimental data set as waveform data using the abfload %function for n = 1:size(directory, 1); x = abfload1(directory(n).name); x0 = []; x1 = []; marks = []; %create 3D array of translated waveform data [m n l] %m=# of samples, n=dimension of samples, l=# of sweeps hfcount = 0; for k = 1:size(x,3); x0(:,1,k) = x(1:700, 1, k); x1(:,1,k) = x0(:,1,k) - mean(x0(1:20,1,k)); x2 = x1(:,:,k); %enter a nested function that will iterate through an array of %data and position the three markers appropriately if k == 1 [M1, M2, M3] = markers(x2); marks(:,1,k) = M1; marks(:,2,k) = M2; marks(:,3,k) = M3; elseif var(x2(1:200)) <= 0.01 marks(:,1,k) = marks(:,1,k-1); marks(:,2,k) = marks(:,2,k-1); marks(:,3,k) = marks(:,3,k-1); if k < 46 %limit high frequency data point count %because points are set to zero hfcount = hfcount+1; hf(hfcount) = k; end else [M1, M2, M3] = markers(x2); marks(:,1,k) = M1; marks(:,2,k) = M2; marks(:,3,k) = M3; end end %create time vector for plotting sweeps (only 1st 7ms) time = linspace(0,7,700)'; %user input analysis options correct = 0; while correct ~= 1 button = questdlg({['Would you like to be able to edit the ',... 'placement of popspike markers for '];... ['Mouse #', strtok((directory(n).name), '.')]},... 'Analysis Option', 'Yes', 'No', 'No'); button = lower(button); if strcmp(button, 'yes') || strcmp(button, 'no') correct = 1; else correct =0; end end %------------------------------------------------------------------- %% Display Population Spike Automation %------------------------------------------------------------------- %plot the markers for each mouse at each sweep if strcmp(button, 'yes') if isempty(button1) button1 = questdlg({'Watch placement of Markers';... ['Keep in mind or record the waveform #(s) ',... 'that appear incorrect'];... ['You will later have the opportunity to ', ... 'modify the waveform(s) that need edits']},... 'Reminder', 'Ok', 'Ok'); end hh = waitbar(0,'Initializing Current progress waitbar...',... 'position',[400 200 275 50]); for i = 1:size(x1,3) pic = figure(1); wave = x1(:,:,i); plot(time',wave); hold on; plot(time(marks(2,1,i)), linspace(-1,2), 'r',... time(marks(2,2,i)), linspace(-1,2),... 'r', time(marks(2,3,i)), linspace(-1,2), 'r'); title({['Mouse #', strtok((directory(n).name), '.')];... ['Waveform #', int2str(i)]}); xlabel('time(ms)'); ylabel('Volts(mV)'); hold off; pause(.2); %update current progress waitbar perc = i/size(x1,3); waitbar(perc,hh,sprintf('%.0f%% Current Progress...',perc*100)) end close(pic); close(hh); %----------------------------------------------------------------- %% Prompt User to Edit Automated PS %----------------------------------------------------------------- button2 = questdlg({['Do you need to edit the placement ',... 'of markers on one or more waveform(s) for:']; ['Mouse #', strtok((directory(n).name), '.'), ' ?']},... 'Analysis Option', 'No', 'Yes', 'No'); button2 = lower(button2); button5 = []; while ~strcmp(button5,'no') %keep a loop while user still needs %to edit the placement of waveform markers %get the range of waveforms to edit markers if strcmp(button2, 'yes') rangetest = 0; while rangetest == 0 prompt = ({['Enter comma separated Range of ',... 'waveform #(s) that need editing Example: 12,19']}); dlg_title = 'Enter comma separated Range of Waveform #(s)'; num_lines = 1; options.Resize='on'; options.WindowStyle='normal'; options.Interpreter='tex'; def = {'1,1'}; Range = inputdlg(prompt,dlg_title,num_lines,def,options); range = Range{1}; range = str2num(range); if length(range) < 2 || length(range) > 2 rangeh = warndlg({'Enter comma separated range';... 'For Example:'; '12,19'}); uiwait(rangeh); continue elseif range(1) == 0 || range(2) > size(x1,3) rangeh = warndlg(['Range should be from 1 to ',... num2str(size(x1,3))]); uiwait(rangeh); continue elseif range(1) > range(2) rangeh = warndlg('For Range [a,b], b must be > a'); uiwait(rangeh); else rangetest = 1; end end end %iterate through the range of waveforms we need to edit if strcmp(button2, 'yes') %initialize waitbar for current progress hh = waitbar(0,'Initializing Current progress waitbar...',... 'position',[400 200 275 50]); for i = range(1):range(2) pic = figure(1); wave = x1(:,:,i); plot(1:length(wave),wave); hold on; plot(marks(2,1,i), linspace(-1,2), 'b',... marks(2,2,i), linspace(-1,2),... 'b', marks(2,3,i), linspace(-1,2), 'b'); title({['Mouse #', strtok((directory(n).name), '.')];... ['Waveform #', int2str(i)]}); xlabel('sample#'); ylabel('Volts(mV)'); button3 = myquestdlg(['Are the markers for the ',... 'current waveform correct?'], 'Analysis Option',... 'Yes', 'No','Quit','Yes'); button3 = lower(button3); if strcmp(button3, 'quit') hold off; break; end if strcmp(button3, 'no') yx = []; if i == range(1) && isempty(button4); button4 = questdlg({'*****';... ['You are going to manually place ',... 'all three markers'];... '1...2...3'; ['You must place each marker',... ' and the order is important']; ... 'Only the x-axis coordinate is relevant'},... 'Important Information', 'OK', 'OK'); end but = 0; while but ~= 3 [xi,yi,click] = ginput(1); if isempty(click) return; elseif click ~= 1 continue; elseif ceil(xi+1) > length(wave) rangeh = warndlg(['Stay in bounds',... '& reselect preveious marker']); uiwait(rangeh); continue; else but = but+1; end xi = ceil(xi); yi = ceil(yi); plot(xi,yi-1:.1:yi+1,'rx') yx(:,but) = [wave(xi+1);xi]; end hold off; marks(:,1,i) = yx(:,1); marks(:,2,i) = yx(:,2); marks(:,3,i) = yx(:,3); end hold off; close(pic); perc = (i-range(1)+1)/(range(2)-range(1)+1); waitbar(perc,hh,sprintf('%.0f%% Current Progress...'... ,perc*100)) end close(hh); end %ask if the user still needs to edit data button5 = questdlg({['Do you want to continue editing ',... 'the placement of markers for:'];... ['Mouse #', strtok((directory(n).name), '.')]},... 'Analysis Option', 'No', 'Yes', 'No'); button5 = lower(button5); end end %------------------------------------------------------------------- %% Consolidate PS Data Collection %------------------------------------------------------------------- %compute relavent calculations that will be used to generate graphs and %will also be used to summarize data POP = []; if size(x,3) >= hf(length(hf)) + 100 for k = 1:hf(length(hf)) + 100; POP(k) = (1/2*(marks(1,1,k) + marks(1,3,k))) - marks(1,2,k); end else for k = 1:size(x,3); POP(k) = (1/2*(marks(1,1,k) + marks(1,3,k))) - marks(1,2,k); end end for k = 1:length(hf); POP(hf(k)) = 0; %make all high freq. stim. samples zero %POP(hf(k)) = nan; %erase all hih freq. stim. samples end POPave = mean(POP(1:32)); POPN = 100*(POP/POPave); POPNave = mean(POPN(112:121)); %create a standard sample length for fitting exponential Y to data if length(POPN) >= 136; POPN_1 = POPN((hf(length(hf)) + 1):136); else POPN_1 = POPN((hf(length(hf)) + 1):length(POPN)); end % Y = POPN_1; Y = Y - min(Y); %normalize Y fitlimit = length(Y); %determine length of fitY % %create the time vectors corresponding to the sample rate, 3 samples per minute. One for %POPN and one for Y and fitY k = 0:1/3:length(POPN)/3; %time vector if length(k) ~= length(POPN) TIME = k(1:length(POPN)); end if length(k) ~= length(Y) k = k(1:length(Y)); end %------------------------------------------------------------------- %% Plot POPN vs TIME & Yfit vs Time & Summar Page %------------------------------------------------------------------- %create plots of POPN vs TIME |&| Y vs Time with fitY vs Time and %create a figure with a summary of data output %plot of POPN vs TIME figure(fig+n); set(gcf,'name',[strtok((directory(n).name), '.'), ' ',... 'POPN vs time(min)']); plot(TIME, POPN,['b','o']); hold on; plot(TIME, POPN,['k','-']); title([strtok((directory(n).name), '.'), ' ',... 'POPN vs time(min)'], 'fontsize', 18, 'fontweight', 'bold'); xlabel('time(min)'); ylabel('Percent Recovery'); text(ceil(hf(1)*1/3)-.5,.05*diff(ylim)+min(ylim),... '\downarrow high frequency stimulation',... 'fontweight', 'bold', 'horizontalalignment', 'left',... 'backgroundcolor', 'y'); text(40,.95*diff(ylim)+min(ylim),['Fortieth min Recovery = '... ,num2str(POPNave,'%2.0f'),'% \downarrow'], 'fontweight',... 'bold', 'horizontalalignment', 'right', 'backgroundcolor', 'y'); hold off; fig = fig + 1; %save plot saveas(gcf,[strtok((directory(n).name), '.'), '_POPNvstime'],'pdf'); %close plot close gcf; %update progress bar step = step +1; perc = (step)/(4*size(directory,1)); perc = num2str(perc,'%0.2f'); perc = str2double(perc); waitbar(perc,h,sprintf('%.0f%% Global Progress...',perc*100)) fitlimit = length(Y); button4 = []; while strcmp(button4, 'no') ~= 1 [curvefit,goodness,fitx] = myfitY(Y,'yes',20,'%PS',fig+n); %sample rate =20seconds/sample Tao = 1/curvefit.beta; %tau in seconds R = abs(goodness.rsquare); %coefficient of determination dfe = goodness.dfe; %degrees of freedom rmse = goodness.rmse; %root mean square error sse = goodness.sse; %sum of squares due to error %plot of Y & fitY vs Time set(gca,'nextplot','replacechildren'); title([strtok((directory(n).name), '.'), ' ',... 'Y and fitY vs Time(sec)'], 'fontsize', 18, 'fontweight', 'bold'); ylabel({'Percent Recovery';... 'post high frequency stimulation'}, 'fontsize', 14); %ask if user wants to edit limit of fit curve button4 = myquestdlg(['Do you need to edit the ',... 'limit of the curve fit?'],... 'Analysis Option', 'No', 'Yes', 'No'); button4 = lower(button4); %have user define limit of fit curve if strcmp(button4, 'yes') % use ginput to get user interaction. Continue until terminated with right click XCropp = []; %initialize clicked coordinates variable click_count = 0; %counts number of times clicked lastpoint = 1; %tracks left click, right click and middle click != 1 pointsY = axis; %assigns axis dimensions pointsY = pointsY(3:4); %keeps only mySignal dimensions while lastpoint == 1 [xi,yi,lastpoint] = ginput(1); %ginput calls input routine line([xi,xi],pointsY,'color','r') %uses last input and draws vertical line in red click_count = click_count+1; %increment clicks XCropp(:,click_count) = xi; %save click coordinates, use later to reference new cropped graph end % select subregion of signal using user input startx = round(min(XCropp(:))/20); %convert from seconds to indices if startx >= 0.5*length(Y); %keep start if only trimming right half of curve startx = 1; end endx = round(max(XCropp(:))/20); %convert from seconds to indices if endx <= 0.5*length(Y); %keep length if only trimming left half of curve endx = length(Y); end Y = Y(startx:endx); limitsCrop = [startx,endx]; % save the cropped limits for future use hold off; end fitY = curvefit.alpha*(1 - exp(-curvefit.beta *(fitx))); end fig = fig + 1; %save plot saveas(gcf,[strtok((directory(n).name), '.'),... '_YandfitYvsTime'],'pdf'); %close plot close gcf; %update progress bar step = step +1; perc = (step)/(4*size(directory,1)); perc = num2str(perc,'%0.2f'); perc = str2double(perc); waitbar(perc,h,sprintf('%.0f%% Global Progress...',perc*100)) %data output summary figure figure(fig+n); set(gcf,'name',[strtok((directory(n).name), '.'), ' ',... 'Output Data']); axis([0,20,0,20]); text(9.54,10,{[strtok((directory(n).name), '.'),... ' DataOutputSummary'];... ['Fortieth Minute Recovery = ',num2str(POPNave, '%2.2f'),'%'];... ['Tau = ', num2str(Tao,'%2.4f'), ' (sec)']; 'Fit Statistics:';... ['Root mean square error = ',num2str(rmse,'%2.4f')];... ['Error sum of squares = ',num2str((sse),'%2.4f')];... ['Coefficient of determination (R^2) = ', num2str(R,'%2.4f')];... ['Degrees of freedom = ', num2str(dfe,'%2.2f')]},... 'horizontalalignment', 'center','fontweight', 'bold',... 'backgroundcolor', 'w', 'verticalalignment', 'middle',... 'fontsize', 24); axis off; %save plot saveas(gcf,[strtok((directory(n).name), '.'),... '_DataOutputSummary'],'pdf'); %close plot close gcf; %update progress bar step = step +1; perc = (step)/(4*size(directory,1)); perc = num2str(perc,'%0.2f'); perc = str2double(perc); waitbar(perc,h,sprintf('%.0f%% Global Progress...',perc*100)) %------------------------------------------------------------------- %% Pack Data into text file for export %------------------------------------------------------------------- %create string array of important data stArray = struct; stArray(n).name = ['mouse_', strtok((directory(n).name), '.')]; stArray(n).waveform = x1; stArray(n).markers = marks; stArray(n).POPN = POPN; stArray(n).POPNave = POPNave; stArray(n).Y = Y; stArray(n).fitY = fitY; stArray(n).Time = TIME; %reformat string array of data that will be exported to txt file w=size(stArray(n).waveform,3); l=size(stArray(n).Y,2); q=length(stArray(n).POPN); stWrite = struct; stWrite(n).name = stArray(n).name; stWrite(n).waveform = reshape(stArray(n).waveform,[700,w]); stWrite(n).markers = reshape(stArray(n).markers(1,:,:),[w,3]); stWrite(n).POPN = reshape(stArray(n).POPN,[q,1]); stWrite(n).POPNave = stArray(n).POPNave; stWrite(n).Y = reshape(stArray(n).Y,[l,1]); stWrite(n).fitY = reshape(stArray(n).fitY,[l,1]); stWrite(n).Time = reshape(stArray(n).Time,[q,1]); %create matrix with all data to be written to text file A = zeros(700,size(x0,3)+11); A(1:length(stWrite(n).Time),1) = stWrite(n).Time; A(1:length(time),2) = time; A(1:length(stWrite(n).markers),3) = stWrite(n).markers(:,1); A(1:length(stWrite(n).markers),4) = stWrite(n).markers(:,2); A(1:length(stWrite(n).markers),5) = stWrite(n).markers(:,3); A(1:length(stWrite(n).POPN),6) = stWrite(n).POPN; A(1,7) = stWrite(n).POPNave; A(1,8) = Tao; A(1,9) = R; A(1:length(stWrite(n).Y),10) = stWrite(n).Y; A(1:length(stWrite(n).fitY),11) = stWrite(n).fitY; %iterate through the wavefroms and append to output matrix for j = 12:size(x0,3)+11 A(:,j) = stWrite(n).waveform(:,j-11); end %create header text array head = []; head1 = ['Time(min),','time(ms),','Marker1,','Marker2,',... 'Marker3,','POPN,','POPNave,','Tau,','corrl,','Y,','fitY,']; head2 = 'sweep#1,'; for j = 2:size(x0,3) head2 = [head2,'sweep#', num2str(j), ',']; end head = [head1, head2]; %create format array format1 = '%s '; for j = 1:size(x0,3)+10 format1 = [format1, '%s ']; end format1 = [format1, '%s']; %open and name and insert header in file for storing output data fid = fopen([strtok((directory(n).name), '.'),... '_OutputData', '.txt'], 'wt'); fprintf(fid, format1, head); fclose(fid); %append output matrix to file dlmwrite([strtok((directory(n).name), '.'), '_OutputData', '.txt'],... A, '-append', 'delimiter', ',', 'roffset',1); %create text file of all analysis data for %recovery and tau if n==1 fids = fopen('Summary of Output Data.txt', 'wt'); fprintf(fids, '%s %s %s', ['subjectID#,','%Recovery,','Tau']); end fprintf(fids, '\n%s %s %s', [[strtok((directory(n).name), '.'),','],... [num2str(POPNave, '%2.2f'),','], [num2str(Tao,'%2.2f'),',']]); %update progress bar step = step +1; perc = (step)/(4*size(directory,1)); perc = num2str(perc,'%0.2f'); perc = str2double(perc); waitbar(perc,h,sprintf('%.0f%% Global Progress...',perc*100)) if n ~= size(directory,1) button6 = myquestdlg(['Do you want to stop ',... 'analyzing the data?'],... 'Analysis Option', 'No', 'Yes', 'No'); button6 = lower(button6); if strcmp(button6, 'yes') break; %n = size(directory, 1)+1; end end end %close progress bar fclose(fids); close(h);