Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

MatlabMain/VSD_rise_delay.m

Go to file
 
 
Cannot retrieve contributors at this time
282 lines (254 sloc) 12.1 KB
Raw Blame
Open in GitHub Desktop
function files = VSD_rise_delay(files,exp_id,rep_id,p_r,p_c,index,darkframe,interval)
%% visualize delays in VSD signals
for ii = exp_id
for jj = rep_id
if isempty(index)
% valid for ANDOR recordings
index{1} = data_rec(ii,jj).metadata.width;
index{2} = data_rec(ii,jj).metadata.height;
index{3} = data_rec(ii,jj).metadata.NumberImages;
end
step_r = 1:p_r:length(index{1});
step_c = 1:p_c:length(index{2});
roi_id = zeros(p_r * p_c * numel(step_r)*numel(step_c) ,1);
pxl_id = reshape(1:length(index{1})*length(index{2}),length(index{1}),length(index{2})); %reshape(1:256*256,256,256);
counter = 0;
for c = step_c
for r = step_r
counter = counter + 1;
tempid = pxl_id(r: r + p_r-1, c: c + p_c-1);
roi_id((1:p_r*p_c) + p_r*p_c*(counter-1) ) = tempid(:);
end
end
% options for subtracting dark frame
if darkframe
temp = files{ii,jj}.RawData.traces{1,2};
else
temp = files{ii,jj}.RawData.traces{1,1};
DF = mean(temp(:,:,1:100),3);
temp = temp - repmat(int16(DF),1,1,size(temp,3));
end
% get the time series for all pixels in rois
v_record = reshape(temp,size(temp,1)*size(temp,2),size(temp,3));
roi_traces = v_record(roi_id,:);
% determine average pixel fluorescence of each roi over time
n_rois = size(roi_traces,1)/(p_r*p_c);
avg_roi_traces = zeros(size(roi_traces,2),n_rois);
for i = 1:size(roi_traces,1)/(p_r*p_c)
indexes = (1:p_r*p_c) + p_r*p_c*(i-1);
tracetemp = mean(roi_traces(indexes,:),1);
avg_roi_traces(:,i) = tracetemp(:);
end
% here find timestamps
start_frame = 1330;
avg_roi_traces = avg_roi_traces(start_frame:end,:);
% % plot rois fluorescence over time
b = size(avg_roi_traces,2); roitest = randi(b,[1 20]);
figure, plot(avg_roi_traces(:,roitest))
% cols = [rand(10,1),rand(10,1),rand(10,1)]
% figure,
% for i=1:5
% plot(squeeze(mean(roi_traces(300:(300+20*i)-1,:),1)),'-','Color',cols(i,:)), hold on
% end
% % legend('1-10','1-20','1-30','1-40','1-50','1-60','1-70','1-80','1-90','1-100')
% legend('1-20','1-40','1-60','1-80','1-100')
% plot(squeeze(mean(roi_traces(300:(300+20*i)-1,:),1)),'-','Color',cols(i,:),'LineWidth',6)
%
xlabel('frame number, 1000 Hz','FontSize',30,'FontWeight','bold')
ylabel('average smoothed and up-sampled dFoF','FontSize',30,'FontWeight','bold')
fig=gcf;
set(findall(fig,'-property','FontSize'),'FontSize',30,'FontWeight','bold')
set(gcf,'color','white')
title('Peak/Event finder','FontSize',50,'FontWeight','bold')
% median filter (fast implementation compiled in C)
F0 = zeros(size(avg_roi_traces));
W = 1500;
tot_length = size(avg_roi_traces,1);
for i = 1:size(avg_roi_traces,2)
x = avg_roi_traces(:,i);
xic = x(1:floor(W/2));
xfc = x(tot_length - floor(W/2)+1:tot_length);
F0(:,i) = fastmedfilt1d(x, W, xic, xfc);
end
dFoF = (avg_roi_traces - F0) ./ F0;
% invert sign of dFoF traces if activity peaks are negatives
mean_dFoF = mean(dFoF,2);
squared_dFoF = mean_dFoF.^2;
activity_times = find(squared_dFoF > prctile(squared_dFoF,[80]));
sign_peaks = sign(mean(mean(dFoF(activity_times,:))));
dFoF = sign_peaks .* dFoF;
% fit a polynomial to smooth traces (alternatives!)
% interval = 866:size(dFoF,1);
s_dFoF = zeros(size(dFoF));
for i = 1:size(dFoF,2)
trace = dFoF(:,i);
s1 = smooth(trace,25,'sgolay',3);
s_dFoF(:,i) = s1;
end
% s_dFoF2 = zeros(size(dFoF));
% s_dFoF3 = zeros(size(dFoF));
% s2 = fastmedfilt1d(trace,50,trace(1:25),trace(end-24:end));
% s_dFoF2(:,i) = s2;
% x = 1:length(trace);
% x = x';
% y = trace;
% s3 = smooth(trace,0.01,'loess');
% s_dFoF(:,i) = s3;
% select rois to test traces and delays
roi_map = reshape(1:length(step_r)*length(step_c),length(step_r),length(step_c));
roitest = diag(roi_map([20:2:30],[5:5:30]));
% up-sample and interpolate the smoothed traces on finer time scale
f = 100; % Hz, frequency data aquisition (1 frame every 10 ms)
f1 = 1000; % 1 frame every 1 ms
x = [1:size(s_dFoF,1)];
xq = [x(1):f/f1:x(end)];
dFoF_interp = interp1(x,s_dFoF,xq,'spline');
% determine time-frame of events from average dFoF
trace = mean(dFoF_interp,2);
% special case of inverted transients
% trace1 = trace - sum([1:length(trace)].*trace')/sum([1:length(trace)]);
% trace2 = sum([1:length(trace)].*trace')/sum([1:length(trace)]) - trace1;
% [pks,locs] = findpeaks(trace2,'MinPeakDistance',300); % get location of local maxima
[pks,locs] = findpeaks(trace,'MinPeakDistance',300); % get location of local maxima
thresh = prctile(trace,90);
ind = find(pks > thresh);
pp = pks(ind);
loc = locs(ind);
m_pp = median(pp);
indd = find(pp>m_pp);
% keep peaks above median of their distribution (filter out spurious secondary peaks)
pp = pp(indd);
loc = loc(indd);
% determine transitions between separate events
cutoff = median(trace);
transitions = NaN(numel(pp)-1,1);
for j = 1:length(pp)-1
avgtracetemp = trace(loc(j):loc(j+1));
transitions(j) = any(avgtracetemp < cutoff);
end
% determine half-rise delays
peakfinal = pp(find(transitions)+1);
locfinal = loc(find(transitions)+1);
% find first peak for each event (events must be separated by at least 1000
% time points)
% crossing1 = zeros(length(step_r),length(step_c),length(peakfinal));
% crossing2 = zeros(length(step_r),length(step_c),length(peakfinal));
crossing3 = zeros(length(step_r),length(step_c),length(peakfinal));
delay_peak = zeros(length(step_r),length(step_c),length(peakfinal));
% location = NaN(size(dFoF_interp,2),length(peakfinal));
% peaks = NaN(size(dFoF_interp,2),length(peakfinal)); % peaks = cell(size(dFoF_interp,2),2);
% half_height = NaN(size(dFoF_interp,2),length(peakfinal)); % half_height = cell(size(dFoF_interp,2),2);
% H_height_delay = NaN(size(dFoF_interp,2),length(peakfinal));% peaks{1,1} = 'peak values';
% peaks{1,2} = 'locations peaks';
% half_height{1,1} = 'half height';
% half_height{1,2} = 'delays rising phase';
% for i = 1:size(dFoF_interp,2)
% trace = dFoF_interp(:,i);
for e = 1:length(peakfinal)
if locfinal(e) - 629 > 0
window = locfinal(e)-629 : locfinal(e)+150;
% window = flipud(window(:));
else
window = 1 : locfinal(j)+150;
% window = flipud(window(:));
end
trace = dFoF_interp(window,:);
% trace11 = mean(trace,2);
% com = sum([1:400].*trace11(1:400)')/sum([1:400]);
% trace12 = trace - repmat(com,size(trace,1),size(trace,2));
% trace = com - trace12;
% temp2 = trace(window);
% b_temp2 = (trace - repmat(min(trace),size(trace,1),1));
% nn_temp2 = b_temp2 ./ repmat(max(b_temp2),size(b_temp2,1),1);
% cutoff = mean(nn_temp2(1:350,:)) + 3*std(nn_temp2(1:350,:));
% % rois_null = find(cutoff > 0.5);
% temp3 = nn_temp2 - repmat(cutoff,size(nn_temp2,1),1);
% [~,min_crossing] = min(abs(temp3));
%
% crossing1(:,:,e) = reshape(min_crossing,length(step_r),length(step_c));
%
% temp4 = nn_temp2 - 0.5; %repmat(cutoff,size(nn_temp2,1),1);
% [~,min_crossing2] = min(abs(temp4));
%
% crossing2(:,:,e) = reshape(min_crossing2,length(step_r),length(step_c));
%% crossing during rising phase of VSD signal based on descen algorithm that localizes half-the-maximum-response
% delays to peak are
[trace_peak peakpos] = sort(trace(end-199:end,:),'descend');
peaks = median(trace_peak(1:5,:));
temp5 = abs( trace-repmat(peaks./2,size(trace,1),1));
[~,min_crossing3] = min(abs(temp5 - 0.5));
crossing3(:,:,e) = reshape(min_crossing3,length(step_r),length(step_c));
peak_CoM = sum(trace_peak(1:5,:).* peakpos(1:5,:),1) ./ ...
sum(trace_peak(1:5,:),1);
peak_delays = round(peak_CoM + length(window)-201);
%
delay_peak(:,:,e) = reshape(peak_delays,length(step_r),length(step_c));
% rng = range(temp2);
% crossing = max(temp2) - rng/2;
% [~, min_crossing] = min(abs(temp2 - crossing));
% H_height_delay(i,jj) = window(min_crossing);
% H_height(i,jj) = temp2(min_crossing);
% [ptemp ltemp] = max(temp2);
% peaks(i,jj) = ptemp;
% location(i,jj) = window(ltemp);
end
% end
% files{ii,jj}.dFoF.data = dFoF;
% files{ii,jj}.dFoF.smoothed = dFoF;
% files{ii,jj}.dFoF.interpolated = [];
delays = H_height_delay - repmat(min(H_height_delay),size(H_height_delay,1),1);
end
end
%% diagnostics plots
x = [25:27;57:59;28:30]'; y = [21:23;25:27;5:7]';
w = max(x(:,1)) - min(x(:,1)); h = max(y(:,1)) - min(y(:,1));
linear_ind = reshape(1:size(step_r,2)*size(step_c,2),size(step_r,2),size(step_c,2));
pairs = [];
for r = 1:size(x,2)
[p,q] = meshgrid(y(:,r), x(:,r));
if r == 1
pairs = [p(:) q(:)];
else
pairs = cat(1,pairs,[p(:) q(:)]);
end
end
avg_downsampled_image = reshape(mean(avg_roi_traces,1),size(step_r,2),size(step_c,2));
idx = sub2ind(size(avg_downsampled_image), pairs(:,1),pairs(:,2));
bw_roi = zeros(size(avg_downsampled_image));
bw_roi(idx) = 1;
figure,
subplot(121), imagesc(avg_downsampled_image), colorbar, axis square, hold on
for r = 1:size(x,2)
rectangle('Position',[x(1,r) y(1,r) w h], 'LineWidth',2, 'EdgeColor','k');
end
title('spatial averaging 3x5 pixels')
subplot(122), imagesc(avg_downsampled_image.*bw_roi), axis square, colorbar,
figure, hold on
cols = ['b','g','r'];
for r = 1:size(x,2)
indtemp = linear_ind(y(:,r),x(:,r));
indtemp = indtemp(:);
plot(trace(:,indtemp),cols(r)),
m(r) = plot(mean(trace(:,indtemp),2),cols(r),'LineWidth',4);
end
intensity_threshold = 2999;
mask_plotting = reshape(mean(avg_roi_traces,1),45,83);
mask_plotting(find(mask_plotting<= intensity_threshold)) = 0;
mask_plotting = logical(mask_plotting);
figure,
for f = 1:5
image = crossing3(:,:,f).*mask_plotting;
% min_scale = prctile( image( find(image(:)>0) ),[5]);
% max_scale = max( image( find(image(:)>0) ),[95]);
subplot(2,3,f), imagesc(image, [400 630]), colorbar, axis tight
title('half-max delays, rise VSD signal')
end
figure,
for f = 1:5
image = delay_peak(:,:,f).*mask_plotting;
% min_scale = prctile( image( find(image(:)>0) ),[5]);
% max_scale = max( image( find(image(:)>0) ),[95]);
subplot(2,3,f), imagesc(image, [630 780]), colorbar, axis tight
title('delays peak VSD signal')
end
end