compact_workflow_andor.m

%% % compact workflow for reading in data for VSD exp %

%% specify folder name where data of interested are saved
recording_dir = '\\storage.laur.corp.brain.mpg.de\Data_2\CoattiAlex\Experiments_2015\X14_20150521\Andor_20150521';
    cd(recording_dir); content_folder = dir;
get_name = 'sampleE';
    folder_ids = find(cell2mat(cellfun(@(x) x==1, {content_folder.isdir},'UniformOutput',0))); 
    select_folder = find(~cellfun('isempty', strfind({content_folder(folder_ids).name},get_name)));
    cd(content_folder(folder_ids(select_folder)).name); list = ls;
for l = 1:length(list)
    file_flag(l) = ~isempty(strfind(list(l,:),'sif'));
end
filenames = list([find(file_flag)],:);
%% select subset of recordings of interest
for d = 1:length(select_folder)   
    file_id = [1:3];
    magnification = [4.*ones(1,3)];
    led_power = [35 35 25];
    exc_wavelength = [405 435 405];
    dye_concentration_um = [150];
    for i = 1:length(file_id)
        filename = filenames(file_id(i),:);  
%         p_r = 1; p_c = 1; darkframe = 0;
        if ~isempty(select_folder)
            absfilepath = strcat(recording_dir,filesep,get_name(d,:),filesep,filename);
        else
            absfilepath = strcat(recording_dir,filesep,filename{i});
        end
        [newdata, meta] = sifread_ale(absfilepath,'read_mode',1);
        % rotation of movie matrix (counteract spooling)
        data = zeros(size(newdata),'uint16');
        for j = 1:size(newdata,3)
            data(:,:,j) = imrotate(newdata(:,:,j),90);
        end
        data_rec(i).data = data;
        data_rec(i).metadata = meta(1);
        clear newdata meta data
        data_rec(i).metadata.objective = magnification(i);
        data_rec(i).metadata.excitation_wavelength = exc_wavelength(i);
        data_rec(i).metadata.led_power = led_power(i);
        data_rec(i).metadata.dye_uM = dye_concentration_um;
    end
end
clearvars -except data_rec

%% inspect raw data (function)  
% inputs: data_rec, npixels (multiple of 2)
npixels = 10;
for i = 3%1:size(data_rec,2)
    disp(['experiment filename >>  ',data_rec(i).metadata.FileName]);
    data = data_rec(i).data;
    % obtain pixels slab
    avg_image = mean(data_rec(i).data,3); 
    figure, subplot(221), imagesc(avg_image), colorbar, title('image of slab')
    % automatic threshold
    [n,bin] =  hist(avg_image(:));  
    subplot(222), bar(bin,n), title('distribution intensities all pixels')
    [n_sorted,sort_order] = sort(n,'descend'); modes = sort(sort_order(1:2)); 
    threshold_intensity = bin(sort_order(find(sort_order > modes(1) & sort_order < modes(2),1,'last')));
    th_mask = avg_image; th_mask(find(avg_image < threshold_intensity)) = 0;
    bw_mask = logical(th_mask); 
        %   figure, imagesc(bw_mask), colorbar 
    % label object and keep object corresponding to slab
    L = bwlabel(bw_mask);
    objects = regionprops(bw_mask,'Area'); [~,obj_ind] = max([objects(:).Area]);
    slab_mask = bw_mask;  
    slab_mask(find(L ~= obj_ind)) = 0; 
    pixels = find(logical(slab_mask));
    subplot(223), imagesc(slab_mask), colorbar, title('slab mask')
    subplot(224), hist(avg_image(pixels)), title('distribution intensity slab pixels')
    % high, medium, low intensity pixels
    [pixels_intensities, pixel_ind] = sort(avg_image(pixels));
    pixels_low = pixels( pixel_ind(find(pixels_intensities>prctile(pixels_intensities,10),1,'first'):...
        find(pixels_intensities>prctile(pixels_intensities,10),1,'first')+npixels-1) ); 
    pixels_high = pixels( pixel_ind(find(pixels_intensities>prctile(pixels_intensities,95),1,'first'):...
        find(pixels_intensities>prctile(pixels_intensities,95),1,'first')+npixels-1) );  %pixels(pixel_ind(end-npixels+1:end));
    pixels_mode = pixels(pixel_ind(find(pixels_intensities>median(pixels_intensities),1,'first'):...
        find(pixels_intensities>median(pixels_intensities),1,'first')+npixels-1));
    pixels_select = [pixels_low, pixels_mode, pixels_high];
    %  Locate epileptiform events in average fluorescence trace over time
    avg_I = mean(reshape(data,size(data,1)*size(data,2),size(data,3)),1);
    nFrames = length(avg_I);
    window = round(1/data_rec(i).metadata.KineticCycleTime * 15); % 1500 at 100 Hz % median filter (fast implementation compiled in C)
    minSBWnd = 300;
    minSBInterval = 300;
    avgSignal= squeeze(avg_I)'; %squeeze(mean(mean(VSData,2),3));
    xic = avgSignal(1:floor(window/2));
    xfc = avgSignal(nFrames - floor(window/2)+1:nFrames);
    avgF0 = squeeze(fastmedfilt1d(avgSignal,window,xic,xfc));
    if ismember(data_rec(i).metadata.excitation_wavelength, [525,550,435,460])
        sign_peaks = -1;
    elseif ismember(data_rec(i).metadata.excitation_wavelength, [490,500,405])
        sign_peaks = +1;
    end
    avgdFoF= sign_peaks*(avgSignal-avgF0)./avgF0;
    [BS,BP,BE,BI]=eventLocator(avgdFoF,1:nFrames,[1;1;1;nFrames],minSBWnd,minSBInterval,'stdThresh',10,'medianWindow',2000);
    % compute dF/F
    f = 1/data_rec(i).metadata.KineticCycleTime; % Hz, frequency data aquisition
    f1 = 1000; % 1 frame every 1 ms, upsampling frequency
    DF = zeros([size(data,3),size(pixels_select)]);
    DFoF = zeros([size(data,3),size(pixels_select)]);
%     sDFoF = zeros([size(data,3),size(pixels_select)]);
%     int_dFoF = zeros([length([1:f/f1:size(data,3)]),size(pixels_select)]);
    tot_length = size(data,3);
    data1 = reshape(data,[size(data,1)*size(data,2),size(data,3)]);
    for ii = 1:size(DF,3)
        for jj = 1:size(DF,2)
            x = double(data1(pixels_select(jj,ii),:));
            xic = x(1:floor(window/2));
            xfc = x((tot_length - floor(window/2)+1):tot_length);
            F0temp = fastmedfilt1d(x, window, xic, xfc);
            temp1 = x'-F0temp;
            temp2 = temp1./F0temp;
%             temp3 = smooth(temp2./F0temp,25,'sgolay',3); 
%             xpoints = [1:size(temp3,1)];
%             xq = [xpoints(1):f/f1:xpoints(end)];
%             temp4 = interp1(xpoints,temp3,xq,'spline');
            DF(:,jj,ii) = temp1 *sign_peaks; DFoF(:,jj,ii) = temp2 *sign_peaks; 
%             sDFoF(:,jj,ii)=temp3; int_dFoF(:,jj,ii) = temp4';
        end
    end
    clear data1
    % plots
    for t = 1:length(BP)
        timewindow(:,t) = round(BP(t)-2*f):round(BP(t)+2*f); % frames in time window peak-2 sec:peak+2sec
    end
    traces = DF; 
    amplitudes = squeeze(prctile(traces(timewindow(:),:,:),95,1));
    baseline = setdiff(floor(window/2):size(traces,1)-floor(window/2),timewindow(:));
    noise = squeeze(sqrt(mean(traces(baseline,:,:).^2,1)));
    SNR = amplitudes ./ noise;
    figure, 
    plot(traces(baseline,:,1)), hold on, plot(mean(traces(baseline,:,1),2),'k','LineWidth',2)
    plot(traces(baseline,:,2)+300), plot(mean(traces(baseline,:,2)+300,2),'k','LineWidth',2)
    plot(traces(baseline,:,3)+600), plot(mean(traces(baseline,:,3)+600,2),'k','LineWidth',2), 
    xlabel('frames baseline'),ylabel('F(t)-F0'),title('fluorescence during baseline for increasing intensities, bottom to top')
    % first figure : compare traces at different intensities
    figure, subplot(221), plot(squeeze(traces(timewindow(:),:,1)),'Color',[0.2 0.5 0.9]), hold on
        plot(mean(traces(timewindow(:),:,1),2),'b','LineWidth',4), 
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('pixel with low intensity','FontSize',18,'FontWeight','bold')
    subplot(222), plot(squeeze(traces(timewindow(:),:,2)),'Color',[0.2 0.9 0.5]), hold on
        plot(mean(traces(timewindow(:),:,2),2),'Color',[0 0.45 0],'LineWidth',4), 
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('pixel with medium intensity','FontSize',18,'FontWeight','bold')
    subplot(223), plot(squeeze(traces(timewindow(:),:,3)),'Color',[0.9 0.4 0.4]); hold on
        plot(mean(traces(timewindow(:),:,3),2),'Color',[0.5 0 0],'LineWidth',4), 
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('pixel with max intensity','FontSize',18,'FontWeight','bold')
    subplot(224), p=plot(squeeze(mean(traces(timewindow(:),:,:),2)),'LineWidth',2), hold on
    set(p(1),'Color','b');set(p(2),'Color',[0 0.45 0]);set(p(3),'Color',[0.5 0 0]);
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('average trace for different intensities','FontSize',18,'FontWeight','bold'),
    legend('low intensity pixels','medium intensity pixels','high intensity pixels')
    figure_handle = gcf;
    allAxesInFigure = findall(figure_handle,'type','axes');
    h_axLabels = get(allAxesInFigure,{'XLabel' 'YLabel'}); 
    set(allAxesInFigure,'FontSize',14,'FontWeight','bold')
    set([h_axLabels{:}],'FontSize',14,'FontWeight','bold')
    % second figure : compare traces with different spatial binning
    traces1 = zeros([size(traces,1),round(size(traces,2)/2),size(traces,3)]);
    for b = 1:round(size(traces,2)/2)
        for bb = 1:size(traces,3)
            traces1(:,b,bb) = mean(traces(:,1:b*2,bb),2);
        end
    end
    figure, subplot(221), p=plot(squeeze(traces1(timewindow(:),:,1))); hold on
    for j = 1:length(p)
        set(p(j),'Color',[0.8 0.8 1]-[0.125 0.125 0]*(j-1))
    end
    legend('2x binning','4x binning','6x binning','8x binning','10x binning');
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('avg pixel with low intensity','FontSize',18,'FontWeight','bold')
    subplot(222), p=plot(squeeze(traces1(timewindow(:),:,2))); hold on
    for j = 1:length(p)
        set(p(j),'Color',[0.55 1 0.7]-[0.125 0.05 0.05]*(j-1))
    end
    legend('2x binning','4x binning','6x binning','8x binning','10x binning');
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('avg pixel with medium intensity','FontSize',18,'FontWeight','bold')
    subplot(223),  p=plot(squeeze(traces1(timewindow(:),:,3))); hold on
    for j = 1:length(p)
        set(p(j),'Color',[1 0.8 0.8]-[0.075 0.125 0.125]*(j-1))
    end
    legend('2x binning','4x binning','6x binning','8x binning','10x binning');
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('avg pixel with max intensity','FontSize',18,'FontWeight','bold')
    subplot(224), 
    % plot SNR versus intensities
    scatter(log10(avg_image(pixels_low(:))),SNR(:,1),'b'), hold on
    plot(mean(log10(avg_image(pixels_low(:)))),mean(SNR(:,1)),'sb','MarkerFaceColor','b','MarkerSize',10)
    scatter(log10(avg_image(pixels_mode(:))),SNR(:,2),[],[0 0.45 0]),
    plot(mean(log10(avg_image(pixels_mode(:)))),mean(SNR(:,2)),'s','Color',[0 0.45 0],'MarkerFaceColor',[0 0.45 0],'MarkerSize',10)
    scatter(log10(avg_image(pixels_high(:))),SNR(:,3),[],'r'),
    plot(mean(log10(avg_image(pixels_high(:)))),mean(SNR(:,3)),'rs','MarkerFaceColor','r','MarkerSize',10)
    xlabel('log10(pixel intensity)','FontSize',18,'FontWeight','bold'), ylabel('SNR','FontSize',18,'FontWeight','bold')
    title('SNR of single pixels with different intensities','FontSize',18,'FontWeight','bold')
    figure_handle = gcf; xLIM = get(gca,'XLim'); set(gca,'xLim',[xLIM(1)-range(xLIM)/2 xLIM(2)+range(xLIM)/2]);
    allAxesInFigure = findall(figure_handle,'type','axes');
    h_axLabels = get(allAxesInFigure,{'XLabel' 'YLabel'}); 
    set(allAxesInFigure,'FontSize',14,'FontWeight','bold')
    set([h_axLabels{:}],'FontSize',14,'FontWeight','bold')
end