https://github.com/simkind/Patch-clamp-analysis
Tip revision: bde5c7399d9f7c789feec0ee26ab5dad4a661d90 authored by simkind on 13 January 2021, 04:40:28 UTC
Create How to use.md
Create How to use.md
Tip revision: bde5c73
APanalysis.m
% abffile = '12d06012.abf'
function [Results] = APanalysis(abffile,graph_on,start,stop,dvdtthreshold,xminmax,channelidx,currentidx)
if nargin ~= 8
error('Not enough input arguments')
end
[~,filename,~] = fileparts(abffile);
[d,si,h]=abfload(abffile);
numsweeps = size(d,3); % number of sweeps
%% For each sweep
% meansweepdata = zeros(numsweeps,size(d,1));
for sweep = 1:numsweeps
Results(sweep).filename = abffile;
data = [];
data = d(:,channelidx,sweep);
timemes = 1:length(data)*(si/1000);
Results(sweep).data = data;
Results(sweep).SI = si;
Results(sweep).currentstimulus = d(:,currentidx,sweep);
%% Get baseline data
[Results(sweep).baseline_potential, Results(sweep).baseline_potentialstd] = baseline(data,start,stop);
Results(sweep).baseline_timerange = [start*(si/1000); stop*(si/1000)]; % need to convert to time.
meantracebl(sweep,1) = Results(sweep).baseline_potential;
Results(sweep).current = mean(d(start:stop,currentidx,sweep));
%% Find peak relative to holding potential
thresh = -10;
[pks, loc] = findpeaks(data,'minpeakheight',thresh,'minpeakdistance',1/(si/1000));
%% CHeck to make sure peaks aren't occuring within 1ms of each other
IPI = diff(loc);%
violates = [];
violates = find(IPI < (1000/si));
if ~isempty(violates)
violates(:,2) = violates + 1;
for i = 1:size(violates,1)
if pks(violates(i,1)) > pks(violates(i,2))
violates(i,1) = 0;
elseif pks(violates(i,1)) < pks(violates(i,2))
violates(i,2) = 0;
elseif pks(violates(i,1)) == pks(violates(i,2))
violates(i,2) = 0;
end
end
violates = violates(:,1) + violates(:,2);
pks(violates) = [];
loc(violates) = [];
end
%% Check if spikes occuring during baseline
withinbl = ismember(loc,start:stop+1/(si/1000));
if length(unique(withinbl)) > 1 % if there are spikes within baseline interval
% calculate new baseline without spikes within it
[Results(sweep).baseline_potential, Results(sweep).baseline_potentialstd, pointstouse]...
= nospikebaseline(data,si,start,stop+1/(si/1000),loc,withinbl);
meantracebl(sweep,1) = Results(sweep).baseline_potential;
Results(sweep).current = mean(d(pointstouse,currentidx,sweep));
% removes peaks that are occuring within bl interval
pks(withinbl) = [];
loc(withinbl) = [];
end
%%
if ~isempty(pks) % if there are spikes
peak_times = loc*(si/1000);
peak_amps = pks-Results(sweep).baseline_potential;
numspikes = length(peak_times);
Results(sweep).num_spikes = numspikes;
Results(sweep).peak_times = peak_times;
Results(sweep).peak_idx = loc;
Results(sweep).peak_amplitudes = pks;
Results(sweep).peak_to_baseline = peak_amps;
% % % % % % % wind = 10/(si/1000); % look into 10ms past peak points
% % % % % % % FastAHP = [];
% % % % % % % for i = 1:numspikes
% % % % % % % pklocation = loc(i);
% % % % % % % if length(data) < (pklocation+wind)
% % % % % % % stoppage = length(data);
% % % % % % % else
% % % % % % % stoppage = pklocation+wind;
% % % % % % % end
% % % % % % % [FastAHP(i,1) FastAHP(i,2)] = min(data(pklocation:stoppage));
% % % % % % % FastAHP(i,2) = FastAHP(i,2) + pklocation;
% % % % % % % end
% % % % % % % Results(sweep).FastAHP_Voltage = FastAHP(:,1);
% % % % % % % Results(sweep).FastAHP_Time = FastAHP(:,2)*(si/1000);
% % % % % % % Results(sweep).FastAHP_Baseline = FastAHP(:,1) - Results(sweep).baseline_potential;
%% Threshold
[Results(sweep).threshold_time, Results(sweep).threshold_amplitude,Results(sweep).threshold_index,...
Results(sweep).dvdtthreshold,Results(sweep).dvdt1,Results(sweep).dvdt2, Results(sweep).dvdt1loc, Results(sweep).dvdt1pks] = ...
SpikeThreshold(data,dvdtthreshold,si,numspikes,xminmax,sweep,0,stop,loc);
% plot the results with the fast AHP and peaks as well
Results(sweep).threshold_baseline = Results(sweep).threshold_amplitude - Results(sweep).baseline_potential;
%%%%%%%%%%%%%%%%%%%%add FastAHP here %%%%%%%%%%%%%%%%%%
[FastAHP] = FastAHPfinder(data,si,Results(sweep).dvdt1loc,dvdtthreshold);
Results(sweep).FastAHP_Voltage = FastAHP(:,1);
Results(sweep).FastAHP_Time = FastAHP(:,2)*(si/1000);
Results(sweep).FastAHP_Baseline = FastAHP(:,1) - Results(sweep).baseline_potential;
if graph_on == 1
fff = figure;
plot(1:length(data),data,'b',loc,pks,'r*')
hold on
plot(1:length(data),data,'b',Results(sweep).threshold_index,Results(sweep).threshold_amplitude,'m*')
xlim([xminmax])
plot(FastAHP(:,2),FastAHP(:,1),'g.')
xlabel(sprintf('Points (multiply by %g for ms)',si/1000 ))
ylabel('Volts')
title(sprintf('Sweep %g File:%s',sweep,filename))
end
%% Spikewidth from Baseline
for i = 1:numspikes
threshold = (peak_amps(i)/2) + Results(sweep).baseline_potential;
pklocation = loc(i);
[Results(sweep).SpikeWidth_Baseline(i,1), ~, ~] = halfwidth(threshold,pklocation,data,si);
end
%% Spikewidth from Threshold
for i = 1:numspikes
pklocation = loc(i); % need to go back to points
if ~isnan(Results(sweep).threshold_amplitude(i))
threshold = ((pks(i) - Results(sweep).threshold_amplitude(i))/2) + Results(sweep).threshold_amplitude(i);
[Results(sweep).SpikeWidth_Threshold(i,1), ~,~] = halfwidth(threshold, pklocation,data,si);
else
Results(sweep).SpikeWidth_Threshold(i,1) = NaN;
end
end
%% Start index for search ADP - if there is spikes, search from fast AHP of that single spike
startsearch = (FastAHP(:,2)+1)*(si/1000); % start search from Fast AHP + 1 bin - this variable is in time, not points
else % if no spikes
numspikes = 0;
%% start index for search ADP - if there is no spike, search ADP from end of baseline interval
startsearch = (stop+1)*(si/1000); % IN TIME, NOT POINTS
%% Plot it if graph_on == 1
if graph_on == 1
fff = figure;
plot(data)
xlabel(sprintf('Points (multiply by %g for ms)',si/1000 ))
ylabel('Volts')
title(sprintf('Sweep %g File:%s',sweep,filename))
end
%% All results based on spikes are blank
Results(sweep).num_spikes = NaN;
Results(sweep).peak_times = NaN;
Results(sweep).peak_idx = NaN;
Results(sweep).peak_amplitudes = NaN;
Results(sweep).peak_to_baseline = NaN;
Results(sweep).FastAHP_Voltage = NaN;
Results(sweep).FastAHP_Time = NaN;
Results(sweep).FastAHP_Baseline = NaN;
Results(sweep).threshold_time = NaN;
Results(sweep).threshold_amplitude = NaN;
Results(sweep).threshold_index = NaN;
Results(sweep).threshold_baseline = NaN;
Results(sweep).dvdtthreshold = NaN;
Results(sweep).dvdt1 = NaN;
Results(sweep).dvdt2= NaN;
Results(sweep).SpikeWidth_Baseline = NaN;
Results(sweep).SpikeWidth_Threshold = NaN;
end
%% Make sure theres only 1 spike
if numspikes > 1
error('Sweep %g contains more than 1 spike',sweep)
end
%% When does ADP return to baseline
if ~isempty(startsearch) %if there is a FAHP, then run next code
[ADPendoutput] = ADPend(data,si,Results(sweep).baseline_potential,startsearch);
Results(sweep).ADPend_time = ADPendoutput.ADPend_time;
Results(sweep).ADPend_amplitude = ADPendoutput.ADPend_amplitude;
Results(sweep).ADPduration = ADPendoutput.ADPduration;
Results(sweep).ADP_AreaUnderCurve = ADPendoutput.ADP_AreaUnderCurve;
Results(sweep).ADPenddata = ADPendoutput;
if ~isempty(ADPendoutput.ADPend_idx) && graph_on == 1
blline = zeros(1,length(data));
blline(:) = Results(sweep).baseline_potential;
xx = zeros(1,length(data));
if ~isempty(pks)
xx(1,FastAHP(:,2)+1) = min(get(gca,'ylim'));
end
xx(1,ADPendoutput.ADPend_idx) = min(get(gca,'ylim')); %-100;
figure(fff)
hold on
plot(blline)
bar(xx,'edgecolor','k','facecolor','k')
xlim([0 length(data)])
end
else
Results(sweep).ADPend_time = [];
Results(sweep).ADPend_amplitude = [];
Results(sweep).ADPduration = [];
Results(sweep).ADP_AreaUnderCurve = [];
Results(sweep).ADPenddata = [];
end
%% Find the Peak of ADP
ADPend_idx = ADPendoutput.ADPend_idx;
if ~isempty(pks) % if there is a spike
FastAHPidx = FastAHP(:,2);
[ADPpeak, ADPloc] = max(data(FastAHPidx+1:ADPend_idx));
ADPloc = ADPloc + FastAHPidx;
if ADPpeak < data(FastAHPidx) % if ADP found is smaller than the fast AHP amplitude, set ADP to fast AHP
ADPloc = FastAHPidx;
ADPpeak = data(FastAHPidx);
end
else % if no spike
[ADPpeak, ADPloc] = max(data(stop+1:ADPend_idx));
ADPloc = ADPloc + stop+1;
end
Results(sweep).ADPpeak_amplitude = ADPpeak;
Results(sweep).ADPpeak_idx = ADPloc;
Results(sweep).ADPpeak_time = ADPloc*(si/1000);
if graph_on == 1
ADPpeakline = zeros(1,length(data));
ADPpeakline(:) = Results(sweep).ADPpeak_amplitude;
xx = zeros(1,length(data));
xx(1,ADPloc) = min(get(gca,'ylim')); %-100;
figure(fff)
hold on
plot(ADPpeakline,'c')
bar(xx,'edgecolor','c','facecolor','c')
xlim([0 length(data)])
end
%% Post-Burst AHP Start
[AHPstartoutput] = PostBurstAHPstart(data,si,Results(sweep).baseline_potential,ADPend_idx);
% Results(sweep).AHPstart_time = AHPstartoutput.AHPstart_time; % the AHPstart time should be the same as the ADP end time
Results(sweep).AHPstart_time = Results(sweep).ADPend_time;
% Results(sweep).AHPpeak_negative_amplitude = AHPstartoutput.AHPpeak_negative_amplitude;
% Results(sweep).AHPpeak_negative_time = AHPstartoutput.AHPpeak_negative_time;
Results(sweep).AHPstartData = AHPstartoutput;
if graph_on == 1
xx = zeros(1,length(data));
xx(1,ADPend_idx) = min(get(gca,'ylim')); %-100;
figure(fff)
hold on
bar(xx,'edgecolor','k','facecolor','k')
xlim([0 length(data)])
end
%% Post-Burst AHP Return
if ~isempty(Results(sweep).AHPstart_time) %if there is AHP start, then run next code
[AHPendoutput] = PostBurstAHPend(data,si,Results(sweep).baseline_potential,Results(sweep).AHPstart_time);
Results(sweep).AHPend_time = AHPendoutput.AHPend_time;
Results(sweep).AHPend_amplitude = AHPendoutput.AHPend_amplitude;
Results(sweep).AHPduration = AHPendoutput.AHPduration;
Results(sweep).AHP_AreaUnderCurve = AHPendoutput.AHP_AreaUnderCurve;
Results(sweep).AHPenddata = AHPendoutput;
if ~isempty(AHPendoutput.AHPend_idx) && graph_on == 1
blline = zeros(1,length(data));
blline(:) = Results(sweep).baseline_potential;
xx = zeros(1,length(data));
xx(1,AHPendoutput.AHPend_idx) = min(get(gca,'ylim')); %-100;
figure(fff)
hold on
plot(blline)
bar(xx,'edgecolor','k','facecolor','k')
xlim([0 length(data)])
end
else
Results(sweep).AHPend_time = [];
Results(sweep).AHPend_amplitude = [];
Results(sweep).AHPduration = [];
Results(sweep).AHP_AreaUnderCurve = [];
Results(sweep).AHPenddata = [];
end
%% Find Peak Negative of AHP
AHPstart_idx = ADPend_idx;
AHPstart_idx = str2double(sprintf('%16.f',AHPstart_idx));
AHPend_idx = AHPendoutput.AHPend_idx;
AHPend_idx = str2double(sprintf('%16.f',AHPend_idx));
[AHPpeak_amplitude, C] = min(data(AHPstart_idx+1:AHPend_idx));
AHPpeak_idx = C+AHPstart_idx;
Results(sweep).AHPpeak_negative_amplitude = AHPpeak_amplitude - Results(sweep).baseline_potential;
Results(sweep).AHPpeak_negative_amplitude_nonbaseline = AHPpeak_amplitude;
Results(sweep).AHPpeak_negative_time = AHPpeak_idx*(si/1000);
Results(sweep).AHPpeak_negative_idx = AHPpeak_idx;
if graph_on == 1
xx = zeros(1,length(data));
xx(1,AHPpeak_idx) = min(get(gca,'ylim')); % -100;
figure(fff)
hold on
bar(xx,'edgecolor','k','facecolor','k')
xlim([0 length(data)])
end
%% Tau
if ~isempty(Results(sweep).AHPpeak_negative_time)
[TAU] = tau(data,si,Results(sweep).baseline_potential,Results(sweep).AHPpeak_negative_time);
Results(sweep).Tau_duration = TAU.tau_duration;
Results(sweep).Tau_time = TAU.tau_time;
Results(sweep).Tau_amplitude = TAU.tau_amplitude;
Results(sweep).Tau_amplitude_baseline = TAU.tau_amplitude_baseline;
if ~isempty(Results(sweep).Tau_time) && graph_on == 1
blline = zeros(1,length(data));
blline(:) = Results(sweep).Tau_amplitude;
xx = zeros(1,length(data));
xx(1,TAU.tau_idx) = min(get(gca,'ylim'));% -100;
figure(fff)
hold on
plot(blline,'r')
bar(xx,'edgecolor','r','facecolor','r')
xlim([0 length(data)])
end
else
Results(sweep).Tau_duration = [];
Results(sweep).Tau_time = [];
Results(sweep).Tau_amplitude = [];
Results(sweep).Tau_amplitude_baseline = [];
end
%% Prep data for mean trace
% meansweepdata(sweep,:) = data - Results(sweep).baseline_potential;
end
% %% Calculate Post-Burst AHP on the mean Trace
% sweep = numsweeps + 1;
% meansweep = mean(meansweepdata,1);
% meansweepbl = mean(meansweep(start:stop));
% if graph_on == 1
% mfff = figure;
% plot(meansweep)
% title('Mean Sweep')
% ylabel('Volts')
% % movegui('northwest')
% xlabel(sprintf('Points (multiply by %g for ms)',si/1000 ))
% end
% Results(1,sweep).data = meansweep;
% Results(1,sweep).baseline_potential = meansweepbl;
% Results(1,sweep).THIS_IS_MEAN_SWEEP = 1;
% Results(1,sweep).SI = si;
% Results(1,sweep).baseline_potentialstd = std(meansweep(start:stop));
% Results(1,sweep).filename = abffile;
% Results(sweep).baseline_timerange = [start*(si/1000); stop*(si/1000)];
% Results(sweep).num_spikes = NaN;
% Results(sweep).peak_times = NaN;
% Results(sweep).peak_idx = NaN;
% Results(sweep).peak_amplitudes = NaN;
% Results(sweep).peak_to_baseline = NaN;
% Results(sweep).FastAHP_Voltage = NaN;
% Results(sweep).FastAHP_Time = NaN;
% Results(sweep).FastAHP_Baseline = NaN;
% Results(sweep).threshold_time = NaN;
% Results(sweep).threshold_amplitude = NaN;
% Results(sweep).threshold_index = NaN;
% Results(sweep).threshold_baseline = NaN;
% Results(sweep).dvdtthreshold = NaN;
% Results(sweep).dvdt1 = NaN;
% Results(sweep).dvdt2= NaN;
% Results(sweep).SpikeWidth_Baseline = NaN;
% Results(sweep).SpikeWidth_Threshold = NaN;
% Results(sweep).SpikeWidth_FirstSpike = NaN;
% Results(sweep).ISI = NaN;
% %% Postburst AHP start for mean sweep
% [mAHPstartoutput] = PostBurstAHPstart(meansweep,si,meansweepbl,stop);
% Results(sweep).AHPstart_time = mAHPstartoutput.AHPstart_time;
% Results(sweep).AHPpeak_negative_amplitude = mAHPstartoutput.AHPpeak_negative_amplitude;
% Results(sweep).AHPpeak_negative_time = mAHPstartoutput.AHPpeak_negative_time;
% Results(sweep).AHP1s_time = mAHPstartoutput.AHP1s_time;
% Results(sweep).AHP1s_mean_amplitude = mAHPstartoutput.AHP1s_mean_amplitude;
% Results(sweep).AHPstartData = mAHPstartoutput;
% if graph_on == 1
% blline = zeros(1,length(data));
% blline(:) = Results(sweep).baseline_potential;
% xx = zeros(1,length(data));
% xx(1,AHPstartoutput.AHPstart_idx) = min(get(gca,'ylim')); %-100;
% xx(1,AHPstartoutput.AHPpeak_negative_idx) = min(get(gca,'ylim')); % -100;
% figure(mfff)
% hold on
% plot(blline)
% bar(xx)
% xlim([0 length(data)])
% end
% %% Postburst AHP Tau for mean sweep
% if ~isempty(Results(sweep).AHPpeak_negative_time)
% [mTAU] = tau(meansweep,si,meansweepbl,Results(sweep).AHPpeak_negative_time);
% Results(sweep).Tau_duration = mTAU.tau_duration;
% Results(sweep).Tau_time = mTAU.tau_time;
% Results(sweep).Tau_amplitude = mTAU.tau_amplitude;
% Results(sweep).Tau_amplitude_baseline = mTAU.tau_amplitude_baseline;
% if ~isempty(Results(sweep).Tau_time) && graph_on == 1
% blline = zeros(1,length(data));
% blline(:) = Results(sweep).Tau_amplitude;
% xx = zeros(1,length(data));
% xx(1,mTAU.tau_idx) = min(get(gca,'ylim'));% -100;
% figure(mfff)
% hold on
% plot(blline,'r')
% bar(xx,'edgecolor','r','facecolor','r')
% xlim([0 length(data)])
% end
% else
% Results(sweep).Tau_duration = [];
% Results(sweep).Tau_time = [];
% Results(sweep).Tau_amplitude = [];
% Results(sweep).Tau_amplitude_baseline = [];
% end
% %% Postburst AHP return for mean sweep
% if ~isempty(Results(sweep).AHPstart_time) %if there is AHP start, then run next code
% [mAHPendoutput] = PostBurstAHPend(meansweep,si,meansweepbl,Results(sweep).AHPstart_time);
% Results(sweep).AHPend_time = mAHPendoutput.AHPend_time;
% Results(sweep).AHPend_amplitude = mAHPendoutput.AHPend_amplitude;
% Results(sweep).AHPduration = mAHPendoutput.AHPduration;
% Results(sweep).AHP_AreaUnderCurve = mAHPendoutput.AHP_AreaUnderCurve;
% Results(sweep).AHPenddata = mAHPendoutput;
% if ~isempty(mAHPendoutput.AHPend_idx) && graph_on == 1
% blline = zeros(1,length(data));
% blline(:) = Results(sweep).baseline_potential;
% xx = zeros(1,length(data));
% xx(1,mAHPendoutput.AHPend_idx) = min(get(gca,'ylim')); %-100;
% figure(mfff)
% hold on
% plot(blline)
% bar(xx)
% xlim([0 length(data)])
% end
% else
% Results(sweep).AHPend_time = [];
% Results(sweep).AHPend_amplitude = [];
% Results(sweep).AHPduration = [];
% Results(sweep).AHP_AreaUnderCurve = [];
% Results(sweep).AHPenddata = [];
% end
end