https://github.com/npmitchell/VisceralOrganMorphogenesisViaCalciumPatternedMuscleConstrictions
Raw File
Tip revision: 3835e0a1faa5853939c2967c867dde5cf49d9887 authored by Noah Mitchell on 25 April 2022, 18:57:35 UTC
added code for figures and dependency code
Tip revision: 3835e0a
figure5_analyze_GCAMP_pulses_posteriorFold.m
% Analyze the GCAMP pulses for width of anterior fold
% NPMitchell 2021


preview = false ;
overwrite = false ;
resDirFn = 'posteriorFoldResults' ;

gutMatlabDir = '/mnt/data/code/gut_matlab/' ;
% datdir = '/Users/npmitchell/Desktop/gut/48YGAL4klarGCAMPGFP/analysis' ;
addpath(fullfile(gutMatlabDir, 'addpath_recurse')) ;
addpath_recurse(gutMatlabDir) 
datdir = '/mnt/data/confocal_data/gut/48YGAL4klarGCAMPGFP/PosteriorFold/analysis/' ;
expts = {'202106221951_e2_1p5x40x_90spf_9sps_12pc_2p5um_hcOilMembrane_1p2usDwell',...
    '202106221951_e3_1p5x40x_90spf_9sps_12pc_2p5um_hcOilMembrane_1p2usDwell', ...
    '202106231410_e1_ch1', ...
    '202106231410_e2_ch1', ...
    '202106231830_e1', ...
    '202106231830_e2', ...
    '202106231830_e3' } ;

% poster frames: 202106211253_e3, 202106211440_e1, 202106211440_e2
% For eLife figure 4 supplement: 202106231830_e1: frame 33 = 16.5 minutes

clipY0s = {[70, 70+282], [70,70+300], ...
    [11,11+251], [27,27+300], ...
    [73,73+283], [24, 24+310], [55, 55+308]} ;
clipYAdjs = {[0, 0], [50,-50]}; 
clipXs = {[1, Inf], [1, Inf], ...
    [1, Inf], [1, Inf], ...
    [1, Inf], [1, Inf], [1, Inf]} ;
% Timestep in MINUTES
dts = [1.5, 1.5, ...
    1.5, 1.5, ...
    1.5, 1.5, 1.5] ;
% Conversion in micron / XY pixels
pix2um = [0.24235971346, 0.24235971346, ...     % 202106221951
    0.24275331665, 0.24275331665, ...           % 202106231410
    0.242359199, 0.242359199, 0.242359199] ;    % 202106231830
 
% TODO: 06211253 --> rotate e1 by +2 deg, e2 by +1 deg in post-flipped frame, redo positions 
% 06211440 --> rotate e2 by +1deg, e1 by -1 deg in post-flipped frame


% DEFINE FOLD TIME AS time of 10 micron indentation on ventral side at
% 16*2.5um = 40 um into the tissue -- this is basically the saggital plane
% fold X position at frame foldT, in pixels
foldXs = [556, 558, ...       % 202106221951 e2, e3
    524,600,...               % 202106211440 e1, e2
    600, 559, 542] ;          % 202106231830 e1, e2, e3
% fold frame # at which position is measured
foldTs = [38, 28, ...         % 202106221951 e2, e3
    18,37,...                 % 202106211440
    22,29,44] ;               % 202106231830
% onset of posterior folding in minutes
t0 = 1.5 * [38, 28, ...         % 202106221951 e2, e3
    18,37,...                 % 202106211440
    22,29,44] ;               % 202106231830

% Fixed xlimits in microns
xlimFix = [-50, 50] ;
ylimFix = [] ;
ylimFixAll = [-20, 30] ;
fixXticks = [-40,-20,0, 20, 40] ;
cLimits = [0, 3] ;

% Filtering of drift deltaX options
capDeltaX = [0, 0, ...
    0, 0, ...
    0, 0, 0] ;
medFiltDeltaXW = [0, 0, ...
    0, 0, ...
    0, 0, 0] ;


% Antp domain is about 40 microns, broadens on dorsal side.

% If drift in X, define it here as [frame#, offsetX] 
xfixed = linspace(-50, 50, 100/0.263) ;
for ee = 1:length(expts)
    fns1 = dir(fullfile(datdir, expts{ee}, 'images', '*c001.png')) ;
    fns2 = dir(fullfile(datdir, expts{ee}, 'images', '*c002.png')) ;
    fns3 = dir(fullfile(datdir, expts{ee}, 'images', '*c003.png')) ;
    
    % For each yrange
    for clipyPairIdx = 1:length(clipYAdjs)
        outfn = sprintf([expts{ee} '_results_Yrange%d.mat'], clipyPairIdx) ;
        if ~exist(fullfile(datdir, outfn), 'file') || overwrite
            if ~exist(fullfile(datdir, outfn), 'file')  
                disp(['results not on disk: ee = ' num2str(ee)])
            else
                disp(['overwriting: ee = ' num2str(ee)])
            end
            % this y range is defined
            clipY = clipY0s{ee} + clipYAdjs{clipyPairIdx} ;
            clipX = clipXs{ee} ;
            foldX = foldXs(ee) ;
            foldT = foldTs(ee) ;

            % Look at differences along y for each timepoint
            ntps = length(fns1) ;
            cmap = viridis(ntps) ;
            timestamps = dts(ee) * (0:ntps-1) ;
            timestamps = timestamps - t0(ee) ;
            
            % Hack
            % if clipyPairIdx == 1
            %     deltaX = zeros(ntps, 1) ;
            % end
            
            clf
            for ii = 1:ntps
                im1 = double(imread(fullfile(fns1(ii).folder, fns1(ii).name))) ;
                im2 = double(imread(fullfile(fns3(ii).folder, fns2(ii).name))) ;
                im3 = double(imread(fullfile(fns3(ii).folder, fns3(ii).name))) ;

                if isfinite(clipX(2))
                    c1 = im1(clipY(1):clipY(2), clipX(1):clipX(2)) ;
                    c2 = im2(clipY(1):clipY(2), clipX(1):clipX(2)) ;
                    c3 = im3(clipY(1):clipY(2), clipX(1):clipX(2)) ;
                else
                    c1 = im1(clipY(1):clipY(2), clipX(1):end) ;
                    c2 = im2(clipY(1):clipY(2), clipX(1):end) ;
                    c3 = im3(clipY(1):clipY(2), clipX(1):end) ;
                end
                
                % mean image
                curr = (c1 + c2 + c3) / 3 ; 
                
                if ii == 50
                    pause(1)
                end
                
                % top-hat TRANSFORM 
                % First do top-hat FILTERING, then subtract
                % -- make an elliptical strel
                d12 = abs(c1 - c2) ; 
                d13 = abs(c1 - c3) ;
                d23 = abs(c2 - c3) ;
                
                % KEEP BIGGER
                % lOpening = 2 ;  % DV extent of opening
                % wOpening = 2 ;  % AP extent of opening
                % se_nlong = round(lOpening / pix2um(ee)) ;
                % se_nwide = round(wOpening / pix2um(ee)) ;
                % % se0 = strel('disk', se_nlong) ;
                % % se_nwide = round(wObjects / pix2um(ee)) ; % 2um diameter / pix2um
                % % se_sampl = round(size(se0.Neighborhood, 2) / se_nwide) ;
                % % seN = se0.Neighborhood(:, 1:se_sampl:end) ;
                % se = strel('arbitrary', ones([se_nlong, se_nwide])) ;
                % hatt1 = imtophat(d12, se) ;
                % d12new = d12 - hatt1 ;
                % figure(1); 
                % subplot(2, 2, 1); imagesc(d12); axis equal; colorbar
                % subplot(2, 2, 2); imagesc(hatt1); axis equal; colorbar
                % subplot(2, 2, 3); imagesc(d12-hatt1); axis equal; colorbar
                % subplot(2, 2, 4); imagesc(hatt1-d12); axis equal; colorbar
                
                % BLUR to KEEP BIGGER
                d12new = imgaussfilt(d12, 0.5 / pix2um(ee)) ;
                d13new = imgaussfilt(d13, 0.5 / pix2um(ee)) ;
                d23new = imgaussfilt(d23, 0.5 / pix2um(ee)) ;
                % d12blur = (d12blur - mean(d12blur) > 0) .* d12blur ;
                
                % KEEP SMALLER
                lOpening = 25 ;  % DV extent of opening
                wOpening = 25 ;  % AP extent of opening
                se_nlong = round(lOpening / pix2um(ee)) ;
                se_nwide = round(wOpening / pix2um(ee)) ;
                % se0 = strel('disk', se_nlong) ;
                % se_nwide = round(wObjects / pix2um(ee)) ; % 2um diameter / pix2um
                % se_sampl = round(size(se0.Neighborhood, 2) / se_nwide) ;
                % seN = se0.Neighborhood(:, 1:se_sampl:end) ;
                se = strel('arbitrary', ones([se_nlong, se_nwide])) ;
                hatt12 = imtophat(d12new, se) ;
                hatt23 = imtophat(d23new, se) ;
                hatt13 = imtophat(d13new, se) ;
                
                
                if preview && mod(ii, 25) == 0
                    figure(1); 
                    subplot(2, 2, 1); imagesc(d12new); axis equal; colorbar
                    subplot(2, 2, 2); imagesc(hatt12); axis equal; colorbar
                    subplot(2, 2, 3); imagesc(d12new-hatt12); axis equal; colorbar
                    subplot(2, 2, 4); imagesc(hatt12-d12new); axis equal; colorbar
                    sgtitle('d12 and hatt12')
                    pause(0.3); 
                    subplot(2, 2, 1); imagesc(d23new); axis equal; colorbar
                    subplot(2, 2, 2); imagesc(hatt23); axis equal; colorbar
                    subplot(2, 2, 3); imagesc(d23new-hatt23); axis equal; colorbar
                    subplot(2, 2, 4); imagesc(hatt23-d23new); axis equal; colorbar
                    sgtitle('d23 and hatt23')
                    pause(0.3); 
                    subplot(2, 2, 1); imagesc(d13new); axis equal; colorbar
                    subplot(2, 2, 2); imagesc(hatt13); axis equal; colorbar
                    subplot(2, 2, 3); imagesc(d13new-hatt13); axis equal; colorbar
                    subplot(2, 2, 4); imagesc(hatt13-d13new); axis equal; colorbar
                    sgtitle('d13 and hatt13')
                    pause(0.3)
                    clf
                    plot(sum(hatt12, 1)); hold on; 
                    plot(sum(hatt23, 1)); hold on; 
                    plot(sum(hatt13, 1)); hold on; 
                    sgtitle('hatt12, hatt23, hatt13')
                    pause(0.3)
                end
                
                if ii == 1 
                    if isfinite(clipX(2))
                        wX = clipX(2)-clipX(1) + 1 ;
                    else
                        wX = size(c1, 2)-clipX(1) + 1 ;
                    end
                    dd = ones(length(fns1), wX) ;
                    ss = ones(length(fns1), wX) ;
                    tophats = ones(length(fns1), wX) ;
                    xshifted = zeros(length(fns1), wX) ;
                    xx = pix2um(ee) * (0:wX-1) ;
                    dinterp = zeros(ntps, size(xfixed, 2)) ; 
                    sinterp = zeros(ntps, size(xfixed, 2)) ;
                    deltaX = zeros(ntps, 1) ;
                elseif clipyPairIdx < 3
                    % Find offset in x wrt previous image
                    assert(any(any(abs(curr-prev))))
                    try
                        [deltaX(ii), ~] = ExtPhaseCorrelation(curr, prev) ;
                        LX = clipXs{ee}(2) - clipXs{ee}(1) ;
                        if abs(deltaX(ii)) > capDeltaX(ee)
                            deltaX(ii) = 0 ;
                        end
                        if deltaX(ii) > LX * 0.5 
                            deltaX(ii) = mod(deltaX(ii), LX) ;
                        elseif deltaX(ii) < -LX* 0.5 
                            deltaX(ii) = -mod(abs(deltaX(ii)), LX) ;                    
                        end
                    catch
                        disp('WARNING: could not compute phase corr!')
                        deltaX(ii) = 0 ;
                    end
                else
                    deltaXfn = sprintf([expts{ee} '_results_Yrange%d.mat'], 2) ;
                    load(fullfile(datdir, deltaXfn), 'deltaX')
                end

                % s1 = sum(c1, 1) ;
                % s2 = sum(c2, 1) ;
                % s3 = sum(c3, 1) ;                 
                % d12 = abs(s1 - s2); 
                % d13 = abs(s1 - s3) ;
                % d23 = abs(s2 - s3) ;
                % d12 = abs(sum(c1 - c2, 1)) ; 
                % d13 = abs(sum(c1 - c3, 1)) ;
                % d23 = abs(sum(c2 - c3, 1)) ;

                % Save as image
                diffDir = sprintf(fullfile(fns1(ii).folder, 'diffs_clipY%02d'), clipyPairIdx) ;
                if ~exist(diffDir, 'dir')
                    mkdir(diffDir)
                end
                fn12 = fullfile(diffDir, [fns1(ii).name(1:end-4) '_d12.png']) ;
                fn23 = fullfile(diffDir, [fns1(ii).name(1:end-4) '_d23.png']) ;
                fn13 = fullfile(diffDir, [fns1(ii).name(1:end-4) '_d13.png']) ;
                assert(max(hatt12(:)) < 255)
                % imwrite(uint16(hatt12 * 2^8), fn12) ;
                % imwrite(uint16(hatt23 * 2^8), fn23) ;
                % imwrite(uint16(hatt13 * 2^8), fn13) ;
                
                % Save transient signal in matrix
                sh12 = sum(hatt12, 1) ;
                sh23 = sum(hatt23, 1) ; 
                sh13 = sum(hatt13, 1) ;
                ddii = mean([sh12; sh23; sh13], 1) ;
                tophats(ii, :) = ddii ;
                
                % Background estimation
                sii = min([sh12; sh23; sh13], [], 1) ;
                ss(ii, :) = sii ;
                dd(ii, :) = ddii - sii ;

                % Save image for phase correlation with next image
                prev = curr ;
            end

            % Median filter deltaX
            if medFiltDeltaXW(ee) > 0
                deltaX = movmedian(deltaX, medFiltDeltaXW(ee)) ;
            end
            
            % Register the position so that (1) stabilized and (2) 0um is anterior fold
            for ii = 1:ntps
                % recall dd and ss for this tp
                ddii = dd(ii, :) ;
                sii = ss(ii, :) ;

                % Save shifted ap positions in matrix
                xshiftii = xx+sum(deltaX(1:ii)) ;
                xshifted(ii, :) = xshiftii ;
            end
            % Now overall offset
            xshiftPix = xshifted(foldT, foldX) ;
            xshifted = xshifted - xshiftPix ;
            
            for ii = 1:ntps
                % recall dd and ss for this tp
                ddii = dd(ii, :) ;
                sii = ss(ii, :) ;

                % Save shifted ap positions in matrix
                xshiftii = xshifted(ii, :) ;
                
                % Interpolate shifted results
                dinterp(ii, :) = interp1(xshiftii, ddii, xfixed) ;
                sinterp(ii, :) = interp1(xshiftii, sii, xfixed) ;

                % Factor in translation of the embryo over the timecourse
                plot(xx+sum(deltaX(1:ii)), dd(ii, :), 'color', cmap(ii, :)) ; hold on;

            end

            % First plot -- raw curves
            set(gcf, 'Visible', 'off')
            xlabel('ap position [$\mu$m]', 'interpreter', 'latex')
            ylabel('GCAMP transient signal [a.u.]', 'interpreter', 'latex')    
            title('Transient GCAMP activity')
            cb = colorbar ;
            caxis([min(timestamps), max(timestamps)])
            ylabel(cb, 'time [min]')
            outfigfn = sprintf(fullfile(datdir, expts{ee}, '00_difference_Yrange%d.png'), clipyPairIdx) ;
            saveas(gcf, outfigfn) 

            % Kymograph -- pre shift
            clf
            imagesc(xx, timestamps, dd)
            xlabel('ap position [$\mu$m]', 'interpreter', 'latex')
            ylabel('time [min]', 'interpreter', 'latex')
            title('Transient GCAMP activity')
            cb = colorbar ;
            ylabel(cb, '$ \delta I $ [a.u.]', 'interpreter', 'latex')
            outfigfn = sprintf(fullfile(datdir, expts{ee}, '01_difference_kymo_Yrange%d.png'), clipyPairIdx) ;
            saveas(gcf, outfigfn)

            % Shifted kymo
            clf
            for row = 1:ntps
                imagesc(xshifted(row, :), timestamps(row), dd(row, :));
                hold on;
            end
            xlabel('stabilized ap position [$\mu$m]', 'interpreter', 'latex')
            ylabel('time [min]', 'interpreter', 'latex')
            title('Transient GCAMP signal', ...
                'interpreter', 'latex')
            xlim(xlimFix)
            if ~isempty(ylimFix)
                ylim(ylimFix)
            else
                ylim([min(timestamps), max(timestamps)])
            end
            cb = colorbar ;
            ylabel(cb, '$ \delta I $ [a.u.]', 'interpreter', 'latex')
            outfigfn = sprintf(fullfile(datdir, expts{ee}, '02_difference_kymoShift_Yrange%d.png'), clipyPairIdx) ;
            saveas(gcf, outfigfn)
            
            
            % PLot deltaX
            clf
            subplot(2, 1, 1)
            plot(timestamps, deltaX, '.-') ;
            ylabel('$\delta x$ [$\mu$m]', 'interpreter', 'latex')
            xlabel('time [min]', 'interpreter', 'latex')
            subplot(2, 1, 2)
            plot(timestamps, cumsum(deltaX), '.-') ;
            ylabel('$\int \delta x$ [$\mu$m]', 'interpreter', 'latex')
            xlabel('time [min]', 'interpreter', 'latex')
            outfigfn = sprintf(fullfile(datdir, expts{ee}, '02b_shift_deltaX_Yrange%d.png'), clipyPairIdx) ;
            saveas(gcf, outfigfn)

            % Plot interpolated results
            clf
            timeGrid = timestamps' .* ones(size(dinterp)) ;
            apGrid = xfixed .* ones(size(dinterp)) ;
            imagesc(xfixed, timestamps, dinterp) ;
            xlim(xlimFix)
            if ~isempty(ylimFix)
                ylim(ylimFix)
            else
                ylim([min(timestamps), max(timestamps)])
            end
            xlabel('stabilized, re-gridded ap position [$\mu$m]', 'interpreter', 'latex')
            ylabel('time [min]', 'interpreter', 'latex')
            title('Transient GCAMP signal [registered frame]', ...
                'interpreter', 'latex')
            cb = colorbar ;
            ylabel(cb, '$ \delta I $ [a.u.]', 'interpreter', 'latex')
            outfigfn = sprintf(fullfile(datdir, expts{ee}, ...
                '03_difference_kymoShiftInterp_Yrange%d.png'), clipyPairIdx) ;
            saveas(gcf, outfigfn)
            xticks(fixXticks)
            ylim(ylimFixAll)
            saveas(gcf, ...
                fullfile(datdir, sprintf('expt%02d_kymo_clipY%d.png', ...
                ee, clipyPairIdx)))

            % Average over time
            % Get t=0 idx
            [~, tidx0] = min(abs(timestamps)) ;
            end5 = find(timestamps - 5 > 0, 1);
            end10 = find(timestamps - 10 > 0, 1);
            end15 = find(timestamps - 15 > 0, 1);
            end20 = find(timestamps - 20 > 0, 1);
            end25 = find(timestamps - 25 > 0, 1);
            end30 = find(timestamps - 30 > 0, 1);
            activity00 = mean(dinterp(1:tidx0, :))' ;
            activity05 = mean(dinterp(tidx0:end5, :))' ;
            activity10 = mean(dinterp(tidx0:end10, :))' ;
            activity15 = mean(dinterp(tidx0:end15, :))' ;
            activity20 = mean(dinterp(tidx0:end20, :))' ;
            activity25 = mean(dinterp(tidx0:end25, :))' ;
            activity30 = mean(dinterp(tidx0:end30, :))' ;
            legends = {'$\langle t<$0 minutes$\rangle$', ...
                '$\langle 0<t<$5 minutes$\rangle$', ...
                '$\langle 0<t<$10 minutes$\rangle$', ...
                '$\langle 0<t<$15 minutes$\rangle$', ...
                '$\langle 0<t<$20 minutes$\rangle$', ...
                '$\langle 0<t<$25 minutes$\rangle$', ...
                '$\langle 0<t<$30 minutes$\rangle$'} ;
            colors = flipud(viridis(7)) ;

            clf
            plot(xfixed, activity00, 'color', colors(1, :)); hold on;
            plot(xfixed, activity05, 'color', colors(2, :))
            plot(xfixed, activity10, 'color', colors(3, :))
            plot(xfixed, activity15, 'color', colors(4, :))
            plot(xfixed, activity20, 'color', colors(5, :))
            plot(xfixed, activity25, 'color', colors(6, :))
            plot(xfixed, activity30, 'color', colors(7, :))
            legend(legends, 'interpreter', 'latex')
            xlabel('stabilized ap position [$\mu$m]', 'interpreter', 'latex')
            ylabel('GCAMP transient activity [a.u.]', 'interpreter', 'latex')    
            title('Transient GCAMP activity', ...
                'interpreter', 'latex')
            xlim(xlimFix)
            outfigfn = sprintf(fullfile(datdir, expts{ee}, '07_difference_meanBgSub_Yrange%d.png'), clipyPairIdx) ;
            saveas(gcf, outfigfn) 

            % Filter in space
            close all
            w1um = round(1 / pix2um(ee)) ;
            plot(xfixed, movmedian(activity00, w1um), 'color', colors(1, :)); 
            hold on;
            plot(xfixed, movmedian(activity05, w1um), 'color', colors(2, :)); 
            plot(xfixed, movmedian(activity10, w1um), 'color', colors(3, :)); 
            plot(xfixed, movmedian(activity15, w1um), 'color', colors(4, :)); 
            plot(xfixed, movmedian(activity20, w1um), 'color', colors(5, :)); 
            plot(xfixed, movmedian(activity25, w1um), 'color', colors(6, :)); 
            plot(xfixed, movmedian(activity30, w1um), 'color', colors(7, :)); 
            legend(legends, 'interpreter', 'latex')
            xlabel('stabilized ap position [$\mu$m]', 'interpreter', 'latex')
            ylabel('GCAMP transient activity [a.u.]', 'interpreter', 'latex')    
            title('Transient GCAMP activity, filtered in space 1 $\mu$m', ...
                'interpreter', 'latex')
            xlim(xlimFix)
            outfigfn = sprintf(fullfile(datdir, expts{ee}, '08_difference_meanBgSubSm_Yrange%d.png'), clipyPairIdx) ;
            saveas(gcf, outfigfn) 
            saveas(gcf, fullfile(datdir, ...
                sprintf('expt%02d_means_clipY%d_%s.png', ...
                ee, clipyPairIdx, expts{ee})))


            % Save result
            save(fullfile(datdir, outfn), 'xx', 'xshifted', 'dd', 'deltaX', ...
                'xfixed', 'dinterp', 'sinterp', 'timeGrid', 'apGrid', ...
                'activity00', 'activity05', 'activity10', ...
                'activity15', 'activity20', 'activity25', 'activity30')
        else
            load(fullfile(datdir, outfn), 'xx', 'xshifted', 'dd', 'deltaX', ...
                'xfixed', 'dinterp', 'sinterp', 'timeGrid', 'apGrid')
            timestamps = unique(timeGrid)' ;
            
            % Average over time
            % Get t=0 idx
            [~, tidx0] = min(abs(timestamps)) ;
            end5 = find(timestamps - 5 > 0, 1);
            end10 = find(timestamps - 10 > 0, 1);
            end15 = find(timestamps - 15 > 0, 1);
            end20 = find(timestamps - 20 > 0, 1);
            end25 = find(timestamps - 25 > 0, 1);
            end30 = find(timestamps - 30 > 0, 1);
            activity00 = mean(dinterp(1:tidx0, :))' ;
            activity05 = mean(dinterp(tidx0:end5, :))' ;
            activity10 = mean(dinterp(tidx0:end10, :))' ;
            activity15 = mean(dinterp(tidx0:end15, :))' ;
            activity20 = mean(dinterp(tidx0:end20, :))' ;
            activity25 = mean(dinterp(tidx0:end25, :))' ;
            activity30 = mean(dinterp(tidx0:end30, :))' ;
            legends = {'$\langle t<$0 minutes$\rangle$', ...
                '$\langle 0<t<$5 minutes$\rangle$', ...
                '$\langle 0<t<$10 minutes$\rangle$', ...
                '$\langle 0<t<$15 minutes$\rangle$', ...
                '$\langle 0<t<$20 minutes$\rangle$', ...
                '$\langle 0<t<$25 minutes$\rangle$', ...
                '$\langle 0<t<$30 minutes$\rangle$'} ;
            colors = flipud(viridis(7)) ;

            clf
            plot(xfixed, activity00, 'color', colors(1, :)); hold on;
            plot(xfixed, activity05, 'color', colors(2, :))
            plot(xfixed, activity10, 'color', colors(3, :))
            plot(xfixed, activity15, 'color', colors(4, :))
            plot(xfixed, activity20, 'color', colors(5, :))
            plot(xfixed, activity25, 'color', colors(6, :))
            plot(xfixed, activity30, 'color', colors(7, :))
            legend(legends, 'interpreter', 'latex')
            xlabel('stabilized ap position [$\mu$m]', 'interpreter', 'latex')
            ylabel('GCAMP transient activity [a.u.]', 'interpreter', 'latex')    
            title('Transient GCAMP activity', ...
                'interpreter', 'latex')
            xlim(xlimFix)
            outfigfn = sprintf(fullfile(datdir, expts{ee}, '07_difference_meanBgSub_Yrange%d.png'), clipyPairIdx) ;
            saveas(gcf, outfigfn) 

            % Filter in space
            close all
            w1um = round(1 / pix2um(ee)) ;
            plot(xfixed, movmedian(activity00, w1um), 'color', colors(1, :)); 
            hold on;
            plot(xfixed, movmedian(activity05, w1um), 'color', colors(2, :)); 
            plot(xfixed, movmedian(activity10, w1um), 'color', colors(3, :)); 
            plot(xfixed, movmedian(activity15, w1um), 'color', colors(4, :)); 
            plot(xfixed, movmedian(activity20, w1um), 'color', colors(5, :)); 
            plot(xfixed, movmedian(activity25, w1um), 'color', colors(6, :)); 
            plot(xfixed, movmedian(activity30, w1um), 'color', colors(7, :)); 
            legend(legends, 'interpreter', 'latex')
            xlabel('stabilized ap position [$\mu$m]', 'interpreter', 'latex')
            ylabel('GCAMP transient activity [a.u.]', 'interpreter', 'latex')    
            title('Transient GCAMP activity, filtered in space 1 $\mu$m', ...
                'interpreter', 'latex')
            xlim(xlimFix)
            outfigfn = sprintf(fullfile(datdir, expts{ee}, '08_difference_meanBgSubSm_Yrange%d.png'), clipyPairIdx) ;
            saveas(gcf, outfigfn) 
            saveas(gcf, fullfile(datdir, ...
                sprintf('expt%02d_means_clipY%d_%s.png', ...
                ee, clipyPairIdx, expts{ee})))

            clf;
            plot(activity30)
            
            % Save result
            save(fullfile(datdir, outfn), 'xx', 'xshifted', 'dd', 'deltaX', ...
                'xfixed', 'dinterp', 'sinterp', 'timeGrid', 'apGrid', ...
                'activity00', 'activity05', 'activity10', ...
                'activity15', 'activity20', 'activity25', 'activity30')
        end
    end
    
end


%% Average all experiments together
expts2include = 1:7;
for clipyPairIdx = 1:2
    colors = define_colors ;
    fixTimeStamps = -21:1.5:31.5 ;  % minutes
    kymoM = zeros(length(xfixed), length(fixTimeStamps)) ;
    nsamples = kymoM ;
    
    % For different averaging
    avgMin = [15, 20, 25] ;
    for avgID = 1:length(avgMin)
        edmy = 1 ;
        disp(['Performing ' num2str(avgMin(avgID)) ' min average'])
        close all
        statAll = zeros(length(xfixed), length(expts2include)) ;
        for ee = expts2include
            disp(['loading ee = ' num2str(ee)])
            outfn = sprintf([expts{ee} '_results_Yrange%d.mat'], clipyPairIdx) ;
            resfn = fullfile(datdir, outfn) ;
            load(resfn, 'xfixed', 'activity15', ...
                'activity20', 'activity25', 'timeGrid', 'dinterp')
            timestamps = unique(timeGrid) ;
            acts = {activity15, activity20, activity25, activity30} ;

            w1um = round(1 / pix2um(ee)) ;
            activity = movmedian(acts{avgID}, w1um) ;
            padd = round(5/pix2um(ee)) ;
            
            % no need for this anymore -- all x are xfixed already
            % keep = xfixed > xlimFix(1) & xfixed < xlimFix(2) ;
            % idx2keep = find(keep) ;
            % normIdx = (idx2keep(1) : idx2keep(1)+padd) ;
            % normIdx = [normIdx, idx2keep(end)-padd:idx2keep] ;
            
            normVal = nanmedian(activity) ;
            bgVal = mean(mink([activity(1:padd); activity(end-padd:end)], 10)) ;
            maxVal = maxk(activity(1+padd:end-padd), 10) ;
            maxVal = mean(maxVal) ;
            anorm = (activity - bgVal) / (maxVal - bgVal); 

            % AVERAGE KYMO
            if avgID == 1
                dnorm = (dinterp - normVal) / (maxVal - normVal) ;
                intrp = griddedInterpolant({timestamps, xfixed'}, dnorm, 'nearest', 'none') ;
                [tt, xx] = ndgrid(fixTimeStamps, xfixed) ;
                newKymo = intrp(tt, xx) ;
                nsamples = nsamples + ~isnan(newKymo)' ;
                newKymo(isnan(newKymo)) = 0 ;
                kymoM = kymoM + newKymo' ;
            end
            
            if avgID == 4
                pause(0.001)
            end
            hold on; 
            plot(xfixed, anorm, 'color', colors(edmy, :)); 
            statAll(:, edmy) = anorm ;
            hold on;

            if ee == expts2include(1)
                xall = xfixed ;
            end
            % min(xfixed(keep))
            % max(xfixed(keep))

            % Collate all results into table
            % keep4all = ismember(xall, xfixed(keep)) ;
            allres(:, edmy) = anorm ;
            edmy = edmy + 1;
        end
        xlabel('ap position from anterior fold [$\mu$m]', 'interpreter', 'latex')
        ylabel('normalized transient GCAMP activity [a.u.]', ...
            'interpreter', 'latex')
        title(sprintf('transient calcium activity, %d min average', avgMin(avgID)), ...
            'interpreter', 'latex')
        resfn = fullfile(datdir, 'gcamp_activity_results') ;
        saveas(gcf, [resfn sprintf('_clipY%d_avg%d.png', clipyPairIdx, avgMin(avgID))])
        saveas(gcf, [resfn sprintf('_clipY%d_avg%d.pdf', clipyPairIdx, avgMin(avgID))])

        
        % Take stats
        close all
        fig = figure('units','centimeters','position',[0,0,6,10]);
        
        % kymo panel
        subplot(2, 1, 1)
        kymoM(~isfinite(kymoM)) = NaN;
        kymoMean = (kymoM - min(kymoM(:))) ./ nsamples  ;
        imagesc(xfixed, fixTimeStamps, kymoMean' ); 
        caxis(cLimits)
        % xlabel('ap position from anterior fold [$\mu$m]', 'interpreter', 'latex')
        ylabel('time from fold onset [min]', 'interpreter', 'latex')
        
        hold on;
        redcolor = colors(2, :) ;
        plot(xlimFix, 0*xlimFix + 5, '--', 'color', redcolor)
        plot(xlimFix, 0*xlimFix - 5, '--', 'color', redcolor)
        plot(xlimFix, 0*xlimFix, '-', 'color', redcolor)
        xticks(fixXticks)
        set(gca, 'xticklabels', [])
        ylim(ylimFixAll)
        colormap(viridis)
            
        % curve panel
        subplot(2, 1, 2)
        avgact = nanmean(statAll, 2) ;
        stdact = nanstd(statAll, 0, 2) ;
        lineProps = {'-','color', colors(1, :)} ;
        factor = 1.0 / max(avgact(:)) ;
        avgact = movmean(avgact, 3) ;
        stdact = movmedian(stdact, 15) ;
        h1=shadedErrorBar(xfixed, avgact*factor, stdact*factor, 'lineProps', lineProps) ;
        % ylim([min(avgact*factor - stdact*factor), max(avgact*factor + stdact*factor)])
        ylim([-0.25, 1.25])
        xticks(fixXticks)
        xlabel('ap position from anterior fold [$\mu$m]', 'interpreter', 'latex')
        ylabel('GCaMP activity [a.u.]', ...
            'interpreter', 'latex')
        title(sprintf('%d min average', avgMin(avgID)), ...
            'interpreter', 'latex')
        resfn = fullfile(datdir, resDirFn, 'gcamp_mean_results') ;
        saveas(gcf, [resfn sprintf('_clipY%d_avg%d.png', clipyPairIdx, avgMin(avgID))])
        saveas(gcf, [resfn sprintf('_clipY%d_avg%d.pdf', clipyPairIdx, avgMin(avgID))])
        
        % Plot mean kymo for first pass
        if avgID == 1
            close all; 
            chelixMap = cubehelix(128,0.43,-0.68,1.3,0.4,[0,1.0],[0,1.0]) ;
            kymoMean = (kymoM - min(kymoM(:))) ./ nsamples  ;
            imagesc(xfixed, fixTimeStamps, kymoMean' ); 
            cb = colorbar;
            caxis(cLimits)
            xlabel('ap position from anterior fold [$\mu$m]', 'interpreter', 'latex')
            ylabel('time from fold onset [min]', 'interpreter', 'latex')
            ylabel(cb,'transient GCAMP activity [a.u.]')
            xticks(fixXticks)
            resfn = fullfile(datdir, resDirFn, 'gcamp_kymo_results') ;
            colormap(viridis)
            saveas(gcf, [resfn sprintf('_clipY%d.png', clipyPairIdx)])
            saveas(gcf, [resfn sprintf('_clipY%d.pdf', clipyPairIdx)])
            colormap(chelixMap)
            saveas(gcf, [resfn sprintf('_clipY%d_cubeHelix.png', clipyPairIdx)])
            saveas(gcf, [resfn sprintf('_clipY%d_cubeHelix.pdf', clipyPairIdx)])
            
        end
    end
end

disp('done')
back to top