% 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 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 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')