Browse Source

Upload matlab scripts to generate paper figures and source data excel files

GRamosT 2 years ago
parent
commit
94e790c07a

+ 173 - 0
figures/MakeCT1NeuronsRFs_Fig5.m

@@ -0,0 +1,173 @@
+close all;
+clear all;
+clc;
+
+%% Define search path.
+parentDir = fileparts(fileparts(mfilename('fullpath')));
+saveFolder = [parentDir filesep 'data'];
+load([saveFolder filesep 'Fig5-CT1-processedTableWithBackgroundOnOff'], 'ProcessedTable');
+
+%%
+figDir = [parentDir filesep 'figures' filesep 'Fig5' filesep];
+if ~exist(figDir, 'dir')
+    mkdir(figDir); 
+end
+%%
+% Stimuli of interest.
+stimSetCellStr = getStimSetCellStr(2);
+%% Split table into cell types.
+tokens = arrayfun(@(x)  regexp(x, '.*CT(\d).*', 'tokens'), ProcessedTable.timeSeriesPath);
+cellTypeNumArray = cellfun(@(x) str2num(x{1}{1}), tokens);
+cellTypeNum = unique(cellTypeNumArray);
+for iNeuron = 1: numel(cellTypeNum)
+    StimTableNeuronCell{iNeuron} = ProcessedTable(cellTypeNumArray == cellTypeNum(iNeuron), :);
+end
+%%
+stimInd = 1;
+nStims = numel(stimSetCellStr);
+for iNeuron = 1: numel(cellTypeNum)
+    for jStim = 1: nStims
+        StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
+    end
+end
+
+for iFly = 1: numel(unique(ProcessedTable.flyInd))
+    flyIndsCell = ProcessedTable.flyInd == iFly;
+    flyStims = ProcessedTable(flyIndsCell, :).stimParamFileName;
+    % Find if stim is within the chosen subset.
+    flyStimInds{iFly} = (cellfun(@(x) find(strcmp(x, stimSetCellStr)), flyStims, 'uni', 0));
+    % If it contains stimuli, check that it includes the 5 chosen.
+    if ~isempty(flyStimInds{iFly})
+        setDiff{iFly} = setdiff(1: 5, [flyStimInds{iFly}{:}]);
+    else
+        setDiff{iFly} = -1;
+    end
+    % For flies with all chosen stimuli order the table accordingly.
+    if isempty(setDiff{iFly})
+        sortFlyIndsTemp = find(flyIndsCell);
+        [~, sortInd] = sort([flyStimInds{iFly}{:}]);
+        sortFlyInds{iFly} = sortFlyIndsTemp(sortInd);
+    end
+end
+
+StimOrderedTable = ProcessedTable(vec([sortFlyInds{:}]), :);
+
+%% Split table into cell types.
+
+tokens = arrayfun(@(x)  regexp(x, '.*CT(\d).*', 'tokens'), StimOrderedTable.timeSeriesPath);
+cellTypeNumArray = cellfun(@(x) str2num(x{1}{1}), tokens);
+cellTypeNum = unique(cellTypeNumArray);
+for iNeuron = 1: numel(cellTypeNum)
+    StimTableNeuronCell{iNeuron} = StimOrderedTable(cellTypeNumArray == cellTypeNum(iNeuron), :);
+end
+%%
+nStims = numel(stimSetCellStr);
+for iNeuron = 1: numel(cellTypeNum)
+    for jStim = 1: nStims
+        StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
+    end
+end
+
+%% Now can use this table cell array into the other functions.
+
+StimTable = ProcessedTable;
+StimSubTablesCell = arrayfun(@(x) getStimSubsetTable(ProcessedTable, ...
+                                                     stimSetCellStr, x), ...
+                             1: numel(stimSetCellStr), 'uni', 0);  
+
+%% Plot tuning curves. Do not discard the on.
+
+alignToOff = true;
+alignToAll = false;
+respQualityThreshold = 0.5;
+rSquareThreshold = 0.2;
+padWithNaNs = 0;
+qualityMeasure = 'responseQuality';
+
+%% Make Figure 5b-d
+cellTypeStr = {'CT1'};
+receptiveFields = plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), -20.0, 0, figDir, [cellTypeStr{:}]);
+%% Export data from Figure 5b-d
+dataFileName = [figDir filesep 'Fig5.xls'];
+sheetNamesCell = {'b-d'};
+
+stimNameStrCell = {'_OFF_horizontal_position_deg', '_OFF_vertical_position_deg', ...
+                   '_ON_horizontal_position_deg', '_ON_vertical_position_deg'};
+               
+for iCellType = 1: numel(cellTypeStr)
+    startRow = 1;
+    sheetName = ['Fig5' sheetNamesCell{iCellType}];
+    for iStim = 1: 4
+        parameterLabel = [cellTypeStr{iCellType} stimNameStrCell{iStim}];
+        DataTable = createRFsTable(receptiveFields{iCellType, iStim}, parameterLabel);
+        writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
+        startRow = startRow + size(DataTable, 1) + 2;
+    end
+end
+%% Inline from plotRFsAlignedWithFitGenotypes
+respQualityThreshold = 0.5;
+rSquareThreshold = 0.2;
+respIndThresh = respQualityThreshold;
+rSqThresh = rSquareThreshold;
+% rSqThresh = -20;
+for iNeuron = 1: size(StimTableNeuronStim(:, 2: end), 1)
+    StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
+    [~, responseIndArrayCell, rSquareArrayCell, fitParamsCell, ...
+      flyIndsCell, ~, ~] = getStimArraysForSameRois(StimTableCell);                                          
+    % Align to fit
+    widthCell{iNeuron} = cellfun(@(x) arrayfun(@(y) 2 * (2 * sqrt(log(2))) * y.c1, x, 'uni', 0), ...
+                        fitParamsCell, 'uni', 1); % The 2 factor is because of bar separation.
+    validIndsCell{iNeuron} = cellfun(@(x, y) x(:) > rSqThresh & y(:) > respIndThresh, ...
+                                     rSquareArrayCell, responseIndArrayCell, 'uni', 0)';
+end
+%%
+validOffInds = cellfun(@(x) x{1} & x{2}, validIndsCell, 'uni', 0);
+validOnInds = cellfun(@(x)  x{3} & x{4}, validIndsCell, 'uni', 0);
+%%
+widthCell = cellfun(@(x) cat(1, x{:}), widthCell, 'uni', 0);
+%%
+xOFFWidth = cellfun(@(x, y) x(1, y), widthCell, validOffInds, 'uni', 0);
+yOFFWidth = cellfun(@(x, y) x(2, y), widthCell, validOffInds, 'uni', 0);
+xONWidth = cellfun(@(x, y) x(3, y), widthCell, validOnInds, 'uni', 0);
+yONWidth = cellfun(@(x, y) x(4, y), widthCell, validOnInds, 'uni', 0);
+allWidthsCell = {xOFFWidth, yOFFWidth, xONWidth, yONWidth};
+%%
+for iW = 1: numel(allWidthsCell)
+    allWidthsCell{iW}(cellfun(@isempty, allWidthsCell{iW}, 'uni', 1)) = {pi};
+end
+%%
+hFig = createPrintFig(15 * [ 1 1/4]);
+colors = brewermap(4, 'Set1');
+colors = colors([1, 2, 4, 3], :);
+for iAx = 1: numel(allWidthsCell)
+    hSubAx(iAx) = subplot(1, 4, iAx);
+    hWidthAx = beanPlot(allWidthsCell{iAx}, cellTypeStr, colors, ...
+                        1, gca, [0-eps 100], 'vertical', 0, 1, 0);
+end
+linkaxes(hSubAx, 'y');
+setFontForThesis(hWidthAx, gcf);
+arrayfun(@prettifyAxes, hWidthAx)
+arrayfun(@offsetAxes, hWidthAx)
+figFileName = ['RF-FWHM-bean'];    
+% figFileName = ['RF-FWHM-bean_rSq-20_respInd0.5'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+
+%% Export data from figure 5bc FWHM distributions
+dataFileName = [figDir filesep 'Fig5.xls'];
+
+stimNameStrCell = {'_OFF_bars_horizontal_deg', '_OFF_bars_vertical_deg', ...
+                   '_ON_bars_horizontal_deg', '_ON_bars_vertical_deg'};
+sheetName = ['Fig5b-d'];
+for iCellType = 1: numel(cellTypeStr)
+    for iStim = 1: 4
+        parameterLabel = [cellTypeStr{iCellType} '_FWHM' stimNameStrCell{iStim}];
+        xlswrite(dataFileName, {parameterLabel}, sheetName, ['A' num2str(startRow)])
+        if allWidthsCell{iStim}{iCellType} == pi % pi was placeholder for no data
+            xlswrite(dataFileName, {''}, sheetName, ['B' num2str(startRow)])
+        else
+            xlswrite(dataFileName, allWidthsCell{iStim}{iCellType}, sheetName, ['B' num2str(startRow)])
+        end
+        startRow = startRow + 1;
+    end
+    startRow = startRow + 1;
+end

+ 269 - 0
figures/MakeFigBarRFs_Fig2bc.m

@@ -0,0 +1,269 @@
+close all;
+clear all;
+clc;
+
+%% Define search path.
+parentDir = fileparts(fileparts(mfilename('fullpath')));
+saveFolder = [parentDir filesep 'data'];
+load([saveFolder filesep 'Fig2-Tm9-processedTableWithBackgroundOnOff'], 'ProcessedTable');
+
+%%
+figDir = [parentDir filesep 'figures' filesep 'Fig2' filesep];
+if ~exist(figDir, 'dir')
+    mkdir(figDir); 
+end
+
+%%
+% Stimuli of interest.
+stimSetCellStr = getStimSetCellStr(2);
+%%
+StimTable = ProcessedTable;
+StimSubTablesCell = arrayfun(@(x) getStimSubsetTable(ProcessedTable, ...
+                                                     stimSetCellStr, x), ...
+                             1: numel(stimSetCellStr), 'uni', 0);
+
+%% Plot tuning curves.
+StimTableCell = StimSubTablesCell(2: end);
+% This extracts the unaligned tuning curves (TCs).
+% Retrieve table data in arrays.
+[tcArrayCell, responseIndArrayCell, rSquareArrayCell, fitParamsCell, ...
+ tcExtCellArray, tcArgExtCellArray] = cellfun(@getTuningCurveArrayFromTable, ...
+                                              StimTableCell, 'uni', 0);
+% Align to fit
+posCell = cellfun(@(x) arrayfun(@(y) y.b1, x, 'uni', 0), fitParamsCell, 'uni', 1);
+widthCell = cellfun(@(x) arrayfun(@(y) 2 * (2 * sqrt(log(2))) * y.c1, x, 'uni', 0), ...
+                    fitParamsCell, 'uni', 1); % The 2 factor is because of bar separation.
+
+lagsCell = cellfun(@(x, y) round(x - size(y, 2) / 2), posCell, tcArrayCell, 'uni', 0);
+barRFs = cellfun(@(x, y) alignRFsToAbsMax(x, 0, y'), tcArrayCell, lagsCell, 'uni', 0);
+validIndsCell = cellfun(@(x, y) x(:) > 0.2 & y(:) > 0.5, ...
+                        rSquareArrayCell, responseIndArrayCell, 'uni', 0)';
+%% Blank epoch was already deleted, now delete vertical bars on edges of screen.
+barRFs{2}(:, end-1: end) = [];
+barRFs{4}(:, end-1: end) = [];
+barRFs{2}(:, 1: 2) = [];
+barRFs{4}(:, 1: 2) = [];
+%%
+offX = barRFs{1}(validIndsCell{1}, :);
+offY = barRFs{2}(validIndsCell{2}, :);
+onX = barRFs{3}(validIndsCell{3}, :);
+onY = barRFs{4}(validIndsCell{4}, :);
+
+% normalize to peak and smooth and normalize again.
+maxNorm = @(x, dim) x ./ max(abs(x), [], dim);
+getRow = @(x, nRow) x(nRow, :);
+normRFs = cellfun(@(x) cell2mat(arrayfun(@(y) smooth(getRow(maxNorm(x,2), y))', 1: size(x, 1), 'uni', 0)'), ...
+                  {offX, offY, onX, onY}, 'uni', 0);
+normRFs = cellfun(@(x) maxNorm(x,2), normRFs, 'uni', 0);
+%%
+% Plot lasagna plots, separate figures for off and on, and reapeated panel
+% is to have a colorbar that does not distort existing axis of interest.
+hFig = createPrintFig(15 * [ 1 2/3]);
+hSubAx = plotLasagnaTracesFromCell(gcf, {offX, offY, offY}, brewermap(6, 'Paired'));
+colorbar
+% arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx(1, :))
+arrayfun(@prettifyAxes, hSubAx)
+arrayfun(@offsetAxes, hSubAx)
+% set(hSubAx, 'Visible', 'off')
+figFileName = ['RFcurves-OFF-XY-lasagna'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+%%
+colorbar
+hFig = createPrintFig(15 * [ 1 2/3]);
+hSubAx = plotLasagnaTracesFromCell(gcf, {onX, onY, onY}, brewermap(6, 'Paired'));
+colorbar
+% arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx(1, :))
+arrayfun(@prettifyAxes, hSubAx)
+arrayfun(@offsetAxes, hSubAx)
+% set(hSubAx, 'Visible', 'off')
+figFileName = ['RFcurves-ON-XY-lasagna'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+%% For nat comms editorial policies source data
+dataFileName = [figDir filesep 'Fig2.xls'];
+DataTable = createRFsTable(offX, 'OFF_horizontal_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig2b')
+DataTable = createRFsTable(offY, 'OFF_vertical_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig2b', 'Range', 'A28:IV100')
+DataTable = createRFsTable(onX, 'ON_horizontal_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig2b', 'Range', 'A56:IV100')
+DataTable = createRFsTable(onY, 'ON_vertical_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig2b', 'Range', 'A84:IV100')
+
+%% This produces panel 2
+
+validWidthCell = getCellArrayFromIndexThreshold(widthCell, validIndsCell', 0);
+hFig = createPrintFig(15 * [ 1/3 1/3]);
+hWidthAx = beanPlot(validWidthCell, {'xOFF' 'yOFF' 'xON' 'yON'}, getONOFFXYColors, ...
+                    1, gca, [0-eps 100], 'vertical', 0, 1, 0);
+setFontForThesis(hWidthAx, gcf);
+arrayfun(@prettifyAxes, hWidthAx)
+arrayfun(@offsetAxes, hWidthAx)
+figFileName = ['RF-FWHM-bean'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+%% Write data for fig2c
+xlswrite(dataFileName, {'FHWM-OFF-bars-horizontal'}, 'Fig2c', 'A1')
+xlswrite(dataFileName, validWidthCell{1}, 'Fig2c', 'B1')
+xlswrite(dataFileName, {'FHWM-OFF-bars-vertical'}, 'Fig2c', 'A2')
+xlswrite(dataFileName, validWidthCell{2}, 'Fig2c', 'B2')
+xlswrite(dataFileName, {'FHWM-ON-bars-horizontal'}, 'Fig2c', 'A3')
+xlswrite(dataFileName, validWidthCell{3}, 'Fig2c', 'B3')
+xlswrite(dataFileName, {'FHWM-ON-bars-vertical'}, 'Fig2c', 'A4')
+xlswrite(dataFileName, validWidthCell{4}, 'Fig2c', 'B4')
+
+%% Range and mean of ON RF fwhm
+min([validWidthCell{3}(:); validWidthCell{4}(:)])
+max([validWidthCell{3}(:); validWidthCell{4}(:)])
+mean([validWidthCell{3}(:); validWidthCell{4}(:)])
+std([validWidthCell{3}(:); validWidthCell{4}(:)])
+
+%%
+threshold = 0.5;
+rSquareThreshold = 0.2;
+qualityMeasure = 'responseQuality';
+widthCell = plotParamsThresholdedDistribution(StimTableCell, threshold, qualityMeasure, figDir, rSquareThreshold);
+%% Do stats for panel 2c
+% Not really normal lly distributed
+cellfun(@lillietest, widthCell);
+% Test offX vs onX and offY vs onY
+[pX, meanXDiff, ~] = permutationTest(widthCell{1}, widthCell{3}, 100000);
+[pY, meanYDiff, ~] = permutationTest(widthCell{2}, widthCell{4}, 100000);
+%% Nested permutations
+flyIndsCell = cellfun(@(y) repelem(y.flyInd, cellfun(@numel, y.responseIndex)), ...
+                      StimTableCell, 'uni', 0);
+for iStim = 1: size(StimTableCell, 2)
+    flyIndsCellArray{iStim} = mat2cell(flyIndsCell{iStim}, ...
+                cellfun(@numel, StimTableCell{iStim}.responseIndex));
+    
+    inputCellArray = flyIndsCellArray{iStim};
+    StimTable = StimTableCell{iStim};
+    
+    nonParamTunings = getRoisFeaturePerChannel(cellfun(@(x) x', StimTable.nonParamTuning, 'uni', 0));
+
+    domainCellArray = arrayfun(@(x) x.domain, nonParamTunings, 'uni', 0);
+    domainSize = cellfun(@numel, domainCellArray);
+    outputCellArray{iStim} = cell2mat(inputCellArray(domainSize == mode(domainSize)));
+
+end
+correctedFlyInds = outputCellArray;
+%%
+goodFlyIndsCell = getCellArrayFromIndexThreshold(correctedFlyInds, validIndsCell', 0);
+
+nFlies = cellfun(@(x) numel(unique(x)), goodFlyIndsCell);
+%% Store the fly indices
+xlswrite(dataFileName, {'FlyID-OFF-bars-horizontal'}, 'Fig2c', 'A6')
+xlswrite(dataFileName, goodFlyIndsCell{1}', 'Fig2c', 'B6')
+xlswrite(dataFileName, {'FlyID-OFF-bars-vertical'}, 'Fig2c', 'A7')
+xlswrite(dataFileName, goodFlyIndsCell{2}', 'Fig2c', 'B7')
+xlswrite(dataFileName, {'FlyID-ON-bars-horizontal'}, 'Fig2c', 'A8')
+xlswrite(dataFileName, goodFlyIndsCell{3}', 'Fig2c', 'B8')
+xlswrite(dataFileName, {'FlyID-ON-bars-vertical'}, 'Fig2c', 'A9')
+xlswrite(dataFileName, goodFlyIndsCell{4}', 'Fig2c', 'B9')
+
+%% Still do the nested permutation
+nROIsPerFlyCell = cellfun(@(x) accumarray(x(:), 1), goodFlyIndsCell, 'uni', 0);
+% Delete flies with no rois to analyze.
+for iStim = 1: numel(nROIsPerFlyCell)
+    nROIsPerFlyCell{iStim}(nROIsPerFlyCell{iStim} == 0) = [];
+end% generate continuos numbering.
+reorderedFlyIndsCell = cellfun(@(x) repelem(1: numel(x), x), nROIsPerFlyCell, 'uni', 0);
+%%
+% samples are cells, but we will permute all cells of a fly together.
+% Test only within orientation differences
+sample1 = 1: numel(nROIsPerFlyCell{1});
+sample2 = numel(sample1) + (1: numel(nROIsPerFlyCell{3}));
+
+permutations = 10000;
+sidedness = 'both';
+plotresult = 1;
+
+[obsDiffX, pValX] = nestetPermutationTest(sample1, sample2, ...
+                                        widthCell([1, 3]), ...
+                                        reorderedFlyIndsCell([1, 3]), ...
+                                        permutations, sidedness, plotresult);
+
+sample1 = 1: numel(nROIsPerFlyCell{2});
+sample2 = numel(sample1) + (1: numel(nROIsPerFlyCell{4}));
+
+[obsDiffY, pValY] = nestetPermutationTest(sample1, sample2, ...
+                                        widthCell([2, 4]), ...
+                                        reorderedFlyIndsCell([2, 4]), ...
+                                        permutations, sidedness, plotresult);
+                                    
+%% Plot raw traces for supplementary figure.
+hFig = createPrintFig(15 * [1 1]);
+
+drawRandomSubset = 1;
+for iStim = 1: size(StimTableCell, 2)
+    
+    inputCellArray = StimTableCell{iStim}.paddedEpochStacks;
+    StimTable = StimTableCell{iStim};
+    % To remove flies with different number of epochs than the majority.
+    nonParamTunings = getRoisFeaturePerChannel(cellfun(@(x) x', StimTable.nonParamTuning, 'uni', 0));
+    
+    domainCellArray = arrayfun(@(x) x.domain, nonParamTunings, 'uni', 0);
+    domainSize = cellfun(@numel, domainCellArray);
+    outputCellArray{iStim} = inputCellArray(domainSize == mode(domainSize));
+
+    allRFTraces = cell2mat(cellfun(@(x) cat(1, x{:})', ...
+                                   outputCellArray{iStim}, 'uni', 0));
+
+    size(allRFTraces)
+    rng('default')
+    % imagesc(allRFTraces(validIndsCell{1}, :))
+    hTracesAx(iStim) = subplot(2, 2, iStim);
+    cla
+    plotInds = find(validIndsCell{iStim});
+    if drawRandomSubset
+        plotInds = plotInds(randperm(10));
+    end
+    barResponseTraces{iStim} = allRFTraces(plotInds, :);
+    timePoints{iStim} = (1: numel(allRFTraces(plotInds(1), :))) / 10; % Data interpolated at 10Hz
+    plot(timePoints{iStim}, bsxfun(@plus, allRFTraces(plotInds, :)', ...
+                0.9 * (1: numel(plotInds))), 'LineWidth', 1, 'Color', 0.25 * [1 1 1]);
+    hold on;
+    axis tight
+    hLines = arrayfun(@(x) refline(0, x), 0.9 * (1: numel(plotInds)));
+    [hLines.Color] = deal([1 1 1] * 0.5);
+    [hLines.LineWidth] = deal(0.5);
+    timePointsPerEpoch = size(allRFTraces, 2) / mode(domainSize);
+    hLines = arrayfun(@(x) line([x, x], ylim), ...
+                      timePointsPerEpoch / 10 * (1/3: 1: mode(domainSize)) + 1/3);
+    [hLines.Color] = deal([1 1 1] * 0.75);
+    [hLines.LineWidth] = deal(0.5);
+    [hLines.LineStyle] = deal(':');
+
+    % Send the reference lines to the back.
+    chH = get(gca,'Children');
+    set(gca,'Children',flipud(chH));
+end
+
+xlabel(hTracesAx(3), 'Time (s)')
+xlabel(hTracesAx(4), 'Time (s)')
+ylabel(hTracesAx(1), '\DeltaF/F_0')
+ylabel(hTracesAx(3), '\DeltaF/F_0')
+
+setFontForThesis(hTracesAx, gcf);
+arrayfun(@prettifyAxes, hTracesAx)
+arrayfun(@offsetAxes, hTracesAx)
+figFileName = ['Fig2-suppl-meanTraces-sortedPositions'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+%% Export data for Supplementary Fig 2
+dataFileName = [figDir filesep 'Fig2.xls'];
+
+sheetName = 'SupFig2';
+stimNameStrCell = {'Tm9_OFF_horizontal_bar_time_seconds', 'Tm9_OFF_vertical_bar_time_seconds', ...
+                   'Tm9_ON_horizontal_bar_time_seconds', 'Tm9_ON_vertical_bar_time_seconds'};
+startRow = 1;
+
+for iStim = 1: 4
+    DataTable = createRFsTable(barResponseTraces{iStim}, stimNameStrCell{iStim});
+    DataTable(:, 1).(stimNameStrCell{iStim}) = timePoints{iStim}(:);
+    writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
+    startRow = startRow + size(DataTable, 1) + 2;
+end
+
+%% Now pick a value an compare stimuli. And do some statistics.
+plotBootstrapCoeffOfVariationTest(widthCell, figDir, threshold, rSquareThreshold)
+%% Variance plot
+plotBarVariance(widthCell, figDir, threshold, rSquareThreshold)
+

File diff suppressed because it is too large
+ 1103 - 0
figures/MakeFigTmNeuronsCDMTableRevisionDoGFitConstraint_Fig7.m


+ 719 - 0
figures/MakeFig_T5Tm9_SpaceTimeTuning_corr_Fig6.m

@@ -0,0 +1,719 @@
+ccc
+%% Define search path.
+parentDir = fileparts(fileparts(mfilename('fullpath')));
+saveFolder = [parentDir filesep 'data'];
+load([saveFolder filesep 'Fig6-T5Tm9-processedTable'], 'ProcessedTable');
+
+%%
+figDir = [parentDir filesep 'figures' filesep 'Fig6' filesep];
+if ~exist(figDir, 'dir')
+    mkdir(figDir); 
+end
+
+%%
+stimSetCellStr = getStimSetCellStr(4);
+
+StimTable = getStimSubsetTable(ProcessedTable, stimSetCellStr, 7);
+
+tuningMethods = StimTable.nonParamTuning{1}(1).tuningMethods;
+
+fieldNamesToUse = {'spacing', 'temporalFrequency'};
+[ stimSize, sortedFieldNameInds, fieldNames, stimParamsFlat ] = ...
+    StimulusDomain_ND( StimTable.stimParams{1}{1}, fieldNamesToUse );
+
+%%
+spacing = [StimTable.stimParams{1}{1}(2:end).spacing];
+stimtrans = [StimTable.stimParams{1}{1}(2:end).stimtrans];
+temporalFreq = unique([stimtrans.mean] ./ spacing);
+spatialFreq = 1 ./ unique(spacing);
+
+%%
+selectedTuningMethod = '1F_amplitude';
+
+tcRoiCell = cellfun(@(x) ...
+    getTuningCurvePerROI(x(1).tc, tuningMethods, '1F_amplitude'), ...
+    StimTable.nonParamTuning, 'uni', 0);
+tm9TCs = [tcRoiCell{:}];
+%%
+tcRoiCell = cellfun(@(x) ...
+    getTuningCurvePerROI(x(2).tc, tuningMethods, '1F_amplitude'), ...
+    StimTable.nonParamTuning, 'uni', 0);
+t5TCs = [tcRoiCell{:}];
+%%
+
+[xGrid, yGrid] = meshgrid(spatialFreq, temporalFreq);
+
+%% These are the maps which average is shown in figure 6g
+respIndex = cell2mat(StimTable.responseIndex')';
+
+validInds = all(respIndex > 0.5, 2);
+figure
+hAx = createSquarishSubplotGrid(sum(validInds));
+iAx = 1;
+for iRoi = find(validInds)'
+    imagesc(hAx(iAx),  t5TCs{iRoi})
+    iAx = iAx + 1;
+end
+figure
+hAx = createSquarishSubplotGrid(sum(validInds));
+iAx = 1;
+for iRoi = find(validInds)'
+    imagesc(hAx(iAx),  tm9TCs{iRoi})
+    iAx = iAx + 1;
+end
+%% Export maps for figure 6g
+% Columns are increasing temporal frequency, rows are increasing spatial
+% frequency. But transpose the array so it matches the spatial and temporal frequency
+% axes. Then columns are increasing spatial frequency and rows increasing
+% temporal frequency.
+[spatialFreqGrid, temporalFreqGrid] = meshgrid(spatialFreq, temporalFreq); % xGrid-space, yGrid-time
+% Create table to append columns with ROI tunint to it.
+DataTable = table(spatialFreqGrid(:), temporalFreqGrid(:), 'VariableNames', {'spatial_frequency_cycles_per_deg', 'temporal_frequency_Hz'});
+
+for iRoi = find(validInds)'
+    tm9Map = tm9TCs{iRoi}';
+    t5Map = t5TCs{iRoi}';
+    DataTable.(sprintf(['ROI%d_Tm9_F1_amplitude'], iRoi)) = tm9Map(:);
+    DataTable.(sprintf(['ROI%d_T5_F1_amplitude'], iRoi)) = t5Map(:);
+end
+
+dataFileName = [figDir filesep 'Fig6.xls'];
+
+sheetName = 'Fig6g';
+writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', 'A1');
+
+%% Figure 6h
+hFig = createPrintFig(10*[1/3 1]);        
+support = [-1-eps 1+eps];
+plotOrientation = 'vertical'; % Can change to plotHorizontalBool.
+halveBean = false;
+useAlpha = true;
+bwOption = true;
+mapCorrs = cellfun(@(x,y) corr2(x, y), t5TCs(validInds), tm9TCs(validInds));
+beanPlot({mapCorrs},{'Tm9-T5'},'k', 1, gca, support, plotOrientation, halveBean);
+% beanPlot({mapCorrs},{'Tm4-T5'},'k', 1, gca, support, plotOrientation, halveBean);
+print(gcf, [figDir 'tuningMap-corr-bean' '.pdf'], '-dpdf', '-r300');
+%% Export pairwise correlation values ffrom figure 6h
+
+labelStrCell = sprintfc(['ROI%d_T5vsTm9_F1_Amplitude_CorrelationCoefficient'], find(validInds));
+for iRow = 1: numel(labelStrCell)
+    xlswrite(dataFileName, labelStrCell(iRow), 'Fig6h', ['A' num2str(iRow)])
+    xlswrite(dataFileName, {mapCorrs(iRow)}, 'Fig6h', ['B' num2str(iRow)])
+end
+%% Shuffle control for correlation of maps for confidence interval in Fig 6h
+validT5TCs = t5TCs(validInds);
+validTm9TCs = tm9TCs(validInds);
+mapCorrs = cellfun(@(x,y) corr2(x, y), validT5TCs, validTm9TCs);
+nShuffles = 1000000;
+% exactShuffles = nchoosek(2 * sum(validInds), sum(validInds)); %even 14
+% gives 40116600 shuffles, so 1e6 is ok approximation
+shuffleMapCorrs = zeros(nShuffles, 1);
+
+for iShuffle = 1: nShuffles
+    shuffleInds = randperm(numel(validTm9TCs));
+    shuffleMapCorrs(iShuffle) = mean(cellfun(@(x,y) corr2(x, y), ...
+                                             validT5TCs, ...
+                                             validTm9TCs(shuffleInds)));
+end
+pValueCorrMaps = sum(shuffleMapCorrs > mean(mapCorrs)) / nShuffles;
+corrConfInt = arrayfun(@(x) quantile(shuffleMapCorrs, x), [0.05, 0.95]);
+meanShuffleCorr = mean(shuffleMapCorrs);
+
+%%
+% Plot spatial tuning by collapsing over temporal frequencies.
+spatialTuning{1} = cell2mat(cellfun(@(x) mean(x, 2)', t5TCs, 'uni', 0)');
+spatialTuning{2} = cell2mat(cellfun(@(x) mean(x, 2)', tm9TCs, 'uni', 0)');
+% Plot temporal tuning by collapsing over spatial frequencies.
+temporalTuning{1} = cell2mat(cellfun(@(x) mean(x, 1), t5TCs, 'uni', 0)');
+temporalTuning{2} = cell2mat(cellfun(@(x) mean(x, 1), tm9TCs, 'uni', 0)');
+%%
+spatialTuning = cellfun(@(x) x./max(x, [], 2), spatialTuning, 'uni', 0);
+temporalTuning = cellfun(@(x) x./max(x, [], 2), temporalTuning, 'uni', 0);
+validInds = all(respIndex > 0.5, 2); % Response quality threshold
+
+%% From here we get mean maps in figure 6g
+selectedTuningMethod = '1F_amplitude';
+createPrintFig(10*[1 1]);
+plotROISpatioTemporalTuning(spatialFreq, temporalFreq, mean(cat(3, t5TCs{validInds}), 3))
+print(gcf, [figDir 'tuningMap-roiMean-0p5respInd' selectedTuningMethod 'T5' '.pdf'], '-dpdf', '-r300');
+
+createPrintFig(10*[1 1]);
+plotROISpatioTemporalTuning(spatialFreq, temporalFreq, mean(cat(3, tm9TCs{validInds}), 3))
+print(gcf, [figDir 'tuningMap-roiMean-0p5respInd' selectedTuningMethod 'Tm9' '.pdf'], '-dpdf', '-r300');
+
+%%
+ProcessedTable = ProcessedTable(ProcessedTable.typeOfProblem == 0, :);
+
+for iFly = 1: numel(unique(ProcessedTable.flyInd))
+    flyIndsCell = ProcessedTable.flyInd == iFly;
+    flyStims = ProcessedTable(flyIndsCell, :).stimParamFileName;
+    % Find if stim is within the chosen subset.
+    flyStimInds{iFly} = (cellfun(@(x) find(strcmp(x, stimSetCellStr)), flyStims, 'uni', 0));
+    % If it contains stimuli, check that it includes the 5 chosen.
+    if ~isempty(flyStimInds{iFly})
+        setDiff{iFly} = setdiff(1: 7, [flyStimInds{iFly}{:}]);
+    else
+        setDiff{iFly} = -1;
+    end
+    % For flies with all chosen stimuli order the table accordingly.
+    if isempty(setDiff{iFly})
+        sortFlyIndsTemp = find(flyIndsCell);
+        [~, sortInd] = sort([flyStimInds{iFly}{:}]);
+        sortFlyInds{iFly} = sortFlyIndsTemp(sortInd);
+    end
+end
+
+StimOrderedTable = ProcessedTable(cat(1, sortFlyInds{:}), :);
+StimOrderedTable = StimOrderedTable(StimOrderedTable.flyInd~=7,:);
+%% Split table into cell types.
+tokens = arrayfun(@(x)  regexp(x, '.*Tm(\d).*', 'tokens'), StimOrderedTable.timeSeriesPath);
+cellTypeNumArray = cellfun(@(x) str2num(x{1}{1}), tokens);
+cellTypeNum = unique(cellTypeNumArray);
+for iNeuron = 1: numel(cellTypeNum)
+    StimTableNeuronCell{iNeuron} = StimOrderedTable(cellTypeNumArray == cellTypeNum(iNeuron), :);
+end
+%%
+stimInd = 1;
+nStims = numel(stimSetCellStr);
+for iNeuron = 1: numel(cellTypeNum)
+    for jStim = 1: nStims
+        StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
+    end
+end
+
+%% Extract the parameters for the analysis of spatial RFs.
+clear amp pos width rsquare responseIndex
+jBar = 1;
+for iStim = 2:5
+    dummyTable = StimTableNeuronStim{1, iStim};
+    for jChannel = 1: 2
+        responseIndex{jChannel}(jBar, :) = cell2mat(cellfun(@(x, y) x(jChannel, 1: end - 1*0), dummyTable.responseIndex, 'uni', 0)');        
+        paramsCell = cellfun(@(x) x(jChannel, 1: end - 1*0), dummyTable.paramTuning, 'uni', 0);
+        amp{jChannel}(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.a1, x, 'uni', 1), paramsCell, 'uni', 0)');
+        pos{jChannel}(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.b1, x, 'uni', 1), paramsCell, 'uni', 0)');
+        width{jChannel}(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.c1, x, 'uni', 1), paramsCell, 'uni', 0)');
+        rsquare{jChannel}(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.gof.rsquare, x, 'uni', 1), paramsCell, 'uni', 0)');      
+    end
+    jBar = jBar + 1;
+end
+pos = cellfun(@(x) 2 * (x - 1), pos, 'uni', 0);
+width = cellfun(@(x) 2 * (2 * sqrt(log(2))) * x, width, 'uni', 0);
+%%
+flyInds = repelem(1:size(StimTableNeuronStim{1}), ... 
+                  cellfun(@(x) size(x, 2), StimTableNeuronStim{1}.responseIndex));
+responseIndexOff = cat(1, responseIndex{1}(1:2, :), responseIndex{2}(1:2, :));
+rSquareOff = cat(1, rsquare{1}(1:2, :), rsquare{2}(1:2, :));
+
+minOffRSquare = cell2mat(cellfun(@(x) min(x(1:2, :)), rsquare, 'uni', 0)');
+minOnRSquare = cell2mat(cellfun(@(x) min(x(3:4, :)), rsquare, 'uni', 0)');
+thresholdRSquare = 0.25;
+minRSquare = min(minOffRSquare);
+minRespInd = min(responseIndexOff);
+isAboveThreshold = minRSquare > thresholdRSquare & minRespInd > 0.3;
+
+%% Figure 6f and supplementary figure 5a
+responseIndexOff = cat(1, responseIndex{1}(1:2, :), responseIndex{2}(1:2, :));
+responseIndexOn = cat(1, responseIndex{1}(3:4, :), responseIndex{2}(3:4, :));
+
+rSquareOff = cat(1, rsquare{1}(1:2, :), rsquare{2}(1:2, :));
+rSquareOn = cat(1, rsquare{1}(3:4, :), rsquare{2}(3:4, :));
+
+minOffRSquare = cell2mat(cellfun(@(x) min(x(1:2, :)), rsquare, 'uni', 0)');
+minOnRSquare = cell2mat(cellfun(@(x) min(x(3:4, :)), rsquare, 'uni', 0)');
+
+plotOFF = true;
+if plotOFF
+    minRespInd = min(responseIndexOff);
+    thresholdRSquare = 0.25;
+    minRSquare = min(minOffRSquare);
+    isAboveThreshold = minRSquare > thresholdRSquare & minRespInd > 0.3;
+else
+    if 0
+        minRespInd = max(responseIndexOn);
+        thresholdRSquare = 0;
+        minRSquare = max(minOnRSquare);
+        isAboveThreshold = minRSquare > thresholdRSquare & minRespInd > 0.3;
+    else
+        minRespInd = min(responseIndexOff);
+        thresholdRSquare = 0.25;
+        minRSquare = min(minOffRSquare);
+        isAboveThreshold = minRSquare > thresholdRSquare & minRespInd > 0.3;
+    end
+end
+
+clear hSubAx
+hFig = createPrintFig(20 * [1 2/3]);
+set(hFig,'color','w');
+
+colorLims = [-1 1] + 0*thresholdRSquare;
+fitParamsCell = {width , pos};
+for jOrientation = 1: 2
+    for iParam = 1: numel(fitParamsCell)
+        iAx = numel(fitParamsCell) * (jOrientation - 1) + iParam;
+        hSubAx(iAx) = subplot(2, numel(fitParamsCell), iAx);
+        if plotOFF
+            x = fitParamsCell{iParam}{1}(jOrientation, isAboveThreshold);
+            y = fitParamsCell{iParam}{2}(jOrientation, isAboveThreshold);
+        else
+            x = fitParamsCell{iParam}{1}(2 + jOrientation, isAboveThreshold);
+            y = fitParamsCell{iParam}{2}(2 + jOrientation, isAboveThreshold);
+        end
+        hScatter = scatter(x, y, [], [1 1 1]*0.25, 'filled', 'MarkerFaceAlpha', 1);
+
+        minVal = min([hSubAx(iAx).XLim(1) hSubAx(iAx).YLim(1)]);
+        maxVal = min([hSubAx(iAx).XLim(2) hSubAx(iAx).YLim(2)]);
+        lsline
+        line([minVal maxVal], [minVal maxVal], 'Color', 0.5 * [1 1 1]);
+        hChildren = get(hSubAx(iAx), 'Children');
+        set(gca, 'Children', [hChildren(end); hChildren(1: end - 1)]);
+        axis square;
+        [corrCoef, p] = corr(x(:), y(:));
+        legStr = ['n = ' num2str(numel(x)) ', corr = ' num2str(corrCoef, 2) ', p = ' num2str(p, 2)];
+        hLeg = legend(hScatter, legStr, 'Location', 'northoutside');
+        hLeg.Box = 'off';
+    end
+end
+linkaxes(hSubAx([1,3]), 'xy')
+linkaxes(hSubAx([2,4]), 'xy')
+hYLabel(1) = ylabel(hSubAx(1), 'T5 OFF RF FWHM Horizontal Bar (deg)');
+hXLabel(1) = xlabel(hSubAx(1), 'Tm9 OFF RF FWHM Horizontal Bar (deg)');
+hYLabel(2) = ylabel(hSubAx(2), 'T5 OFF RF Center Horizontal Bar (deg)');
+hXLabel(2) = xlabel(hSubAx(2), 'Tm9 OFF RF Center Horizontal Bar (deg)');
+hYLabel(3) = ylabel(hSubAx(3), 'T5 OFF RF FWHM Vertical Bar (deg)');
+hXLabel(3) = xlabel(hSubAx(3), 'Tm9 OFF RF FWHM Vertical Bar (deg)');
+hYLabel(4) = ylabel(hSubAx(4), 'T5 OFF RF Center Vertical Bar (deg)');
+hXLabel(4) = xlabel(hSubAx(4), 'Tm9 OFF RF Center Vertical Bar (deg)');
+
+[hXLabel.FontSize] = deal(10);
+[hYLabel.FontSize] = deal(10);
+hLeg.FontSize = 8;
+[hSubAx.FontSize] = deal(8);
+[hSubAx.Color] = deal('none');
+
+arrayfun(@prettifyAxes, [hSubAx]); 
+if plotOFF
+    print(gcf, [figDir 'rfFitParamsScatter-Off-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');
+else
+    print(gcf, [figDir 'rfFitParamsScatter-On-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');
+end
+%% Export data for figure 6f ans supplementary figure 6a
+tm9XFWHM = width{1}(1, isAboveThreshold);
+t5XFWHM = width{2}(1, isAboveThreshold);
+tm9YFWHM = width{1}(2, isAboveThreshold);
+t5YFWHM = width{2}(2, isAboveThreshold);
+
+tm9XPos = pos{1}(1, isAboveThreshold);
+t5XPos = pos{2}(1, isAboveThreshold);
+tm9YPos = pos{1}(2, isAboveThreshold);
+t5YPos = pos{2}(2, isAboveThreshold);
+
+
+DataTable = table(tm9XFWHM(:), t5XFWHM(:), tm9YFWHM(:), t5YFWHM(:), ...
+                  tm9XPos(:), t5XPos(:), tm9YPos(:), t5YPos(:), ...
+                  'VariableNames', ...
+                   {'Tm9_OFF_RF_FWHM_Horizontal_Bar_deg', ...
+                    'T5_OFF_RF_FWHM_Horizontal_Bar_deg', ...
+                    'Tm9_OFF_RF_FWHM_Vertical_Bar_deg', ...
+                    'T5_OFF_RF_FWHM_Vertical_Bar_deg', ...
+                    'Tm9_OFF_RF_Position_Horizontal_Bar_deg', ...
+                    'T5_OFF_RF_Position_Horizontal_Bar_deg', ...
+                    'Tm9_OFF_RF_Position_Vertical_Bar_deg', ...
+                    'T5_OFF_RF_Position_Vertical_Bar_deg', ...
+                   });
+   
+
+dataFileName = [figDir filesep 'Fig6.xls'];
+
+sheetName = 'Fig6f';
+writetable(DataTable(:, 1: 4), dataFileName, 'Sheet', sheetName, 'Range', 'A1');
+
+sheetName = 'SuppFig6ab';
+writetable(DataTable(:, 5: 8), dataFileName, 'Sheet', sheetName, 'Range', 'A1');
+
+%% Reply to reviwer#1 comment 5
+doBootstrap = 0;
+
+xH = vec(width{1}(1, isAboveThreshold));
+yH = vec(width{2}(1, isAboveThreshold));
+xV = vec(width{1}(2, isAboveThreshold));
+yV = vec(width{2}(2, isAboveThreshold));
+tmpInds = flyInds(isAboveThreshold);
+nValidFlies = numel(unique(tmpInds));
+iAx = 1;
+
+
+clear hSubAx
+mSize = 15;
+hFig = createPrintFig(25 * [1 1/3]);
+for iFly = unique(tmpInds)   
+    hSubAx(iAx, 1) = subplot(2, nValidFlies, iAx);
+    x = xH(iFly == tmpInds);
+    y = yH(iFly == tmpInds);
+    scatter(xH, yH, mSize, 'filled', 'MarkerFaceAlpha', 0.1, 'MarkerFaceColor', 'k');
+    hold on;
+    hScatter = scatter(x, y, mSize);
+    hLines = lsline;
+    delete(hLines(2));
+%     refline(1,0)
+    if doBootstrap && numel(x) > 1
+        [bootstrapCorrCI, bootstrapCorrs] = bootci(10000, @corr, x(:), y(:));
+        legStr = sprintf('n = %d\ncorr = %.2f\n0.95-CI = [%.2f, %.2f]', ...
+                          numel(x), mean(bootstrapCorrs), bootstrapCorrCI);
+    else
+        [corrCoef, p] = corr(x(:), y(:));
+        legStr = ['n = ' num2str(numel(x)) ...
+                  ', \n corr = ' num2str(corrCoef, 2) ', \n p = ' num2str(p, 2)];
+    %     hLeg = legend(hScatter, legStr, 'Location', 'northoutside');
+    end
+%     hLeg.Box = 'off';    
+    title(sprintf(legStr));
+    corrArray(iAx, 1) = corrCoef;
+
+    hSubAx(iAx, 2) = subplot(2, nValidFlies, nValidFlies + iAx);
+    x = xV(iFly == tmpInds);
+    y = yV(iFly == tmpInds);
+    scatter(xV, yV, mSize, 'filled', 'MarkerFaceAlpha', 0.1, 'MarkerFaceColor', 'k');
+    hold on;
+    hScatter = scatter(x, y, mSize);
+    hLines = lsline;
+    delete(hLines(2));
+    if doBootstrap && numel(x) > 1
+        [bootstrapCorrCI, bootstrapCorrs] = bootci(10000, @corr, x(:), y(:));
+        legStr = sprintf('n = %d\ncorr = %.2f\n0.95-CI = [%.2f, %.2f]', ...
+                          numel(x), mean(bootstrapCorrs), bootstrapCorrCI);
+    else
+        [corrCoef, p] = corr(x(:), y(:));
+        legStr = ['n = ' num2str(numel(x)) ...
+                  '\n corr = ' num2str(corrCoef, 2) '\n p = ' num2str(p, 2)];
+    %     hLeg = legend(hScatter, legStr, 'Location', 'northoutside');
+    end
+    title(sprintf(legStr))
+    corrArray(iAx, 2) = corrCoef;
+    iAx = iAx + 1;
+end
+
+axis(hSubAx, 'square');
+linkaxes(hSubAx, 'xy')
+clear hXLabel hYLabel
+hXLabel(1) = xlabel(hSubAx(1, 1), 'Tm9 fwhm hor. (deg)');
+hYLabel(1) = ylabel(hSubAx(1, 1), 'T5 fwhm hor. (deg)');
+hXLabel(2) = xlabel(hSubAx(1, 2), 'Tm9 fwhm vert. (deg)');
+hYLabel(2) = ylabel(hSubAx(1, 2), 'T5 fwhm vert. (deg)');
+
+[hXLabel(1).FontSize] = deal(10);
+[hYLabel(1).FontSize] = deal(10);
+% hLeg.FontSize = 8;
+[hSubAx.FontSize] = deal(8);
+[hSubAx.Color] = deal('none');
+arrayfun(@prettifyAxes, [hSubAx]); 
+
+print(gcf, [figDir 'rfFwhmScatterPerFly-Off-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');
+
+%% average values within fly and plot
+% fly means
+xHFlyMeans = accumarray(tmpInds', xH, [], @mean);
+yHFlyMeans = accumarray(tmpInds', yH, [], @mean);
+xHFlyStds = accumarray(tmpInds', xH, [], @std);
+yHFlyStds = accumarray(tmpInds', yH, [], @std);
+xHFlyMins = accumarray(tmpInds', xH, [], @min);
+yHFlyMins = accumarray(tmpInds', yH, [], @min);
+xHFlyMaxs = accumarray(tmpInds', xH, [], @max);
+yHFlyMaxs = accumarray(tmpInds', yH, [], @max);
+
+[x, y] = deal(xHFlyMeans(xHFlyMeans~=0), yHFlyMeans(xHFlyMeans~=0));
+[xStd, yStd] = deal(xHFlyStds(xHFlyMeans~=0), yHFlyStds(xHFlyMeans~=0));
+[xMin, yMin] = deal(xHFlyMins(xHFlyMeans~=0), yHFlyMins(xHFlyMeans~=0));
+[xMax, yMax] = deal(xHFlyMaxs(xHFlyMeans~=0), yHFlyMaxs(xHFlyMeans~=0));
+
+hFig = createPrintFig(7 * [1 1]);
+hScatter = errorbar(x, y, xMin - x, xMax - x, yMin - y, yMax - y, ...
+                    'LineStyle', 'none', 'CapSize', 0, 'LineWidth', 1, ...
+                    'Color', [1 1 1] * 0.4);
+hold on
+hScatter = scatter(x, y, 25, 'filled', 'MarkerFaceColor', [1 1 1] * 0.20);
+lsline
+if doBootstrap && numel(x) > 1
+    [bootstrapCorrCI, bootstrapCorrs] = bootci(10000, @corr, x(:), y(:));
+    legStr = sprintf('n = %d, corr = %.2f\n0.95-CI = [%.2f, %.2f]', ...
+                      numel(x), mean(bootstrapCorrs), bootstrapCorrCI);
+else
+    [corrCoef, p] = corr(x(:), y(:));
+    legStr = ['n = ' num2str(numel(x)) ...
+              ', \n corr = ' num2str(corrCoef, 2) ', \n p = ' num2str(p, 2)];
+end
+title(sprintf(legStr));
+xlabel('Tm9 fwhm hor. (deg)');
+ylabel('T5 fwhm hor. (deg)');
+% hLeg = legend(hScatter, legStr, 'Location', 'northoutside');
+% hLeg.Box = 'off';   
+prettifyAxes(gca);
+
+print(gcf, [figDir 'rfFwhmScatterFlyMeans-XOff-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');
+%% Same fo vertical
+xHFlyMeans = accumarray(tmpInds', xV, [], @mean);
+yHFlyMeans = accumarray(tmpInds', yV, [], @mean);
+xHFlyStds = accumarray(tmpInds', xV, [], @std);
+yHFlyStds = accumarray(tmpInds', yV, [], @std);
+xHFlyMins = accumarray(tmpInds', xV, [], @min);
+yHFlyMins = accumarray(tmpInds', yV, [], @min);
+xHFlyMaxs = accumarray(tmpInds', xV, [], @max);
+yHFlyMaxs = accumarray(tmpInds', yV, [], @max);
+
+[x, y] = deal(xHFlyMeans(xHFlyMeans~=0), yHFlyMeans(xHFlyMeans~=0));
+[xStd, yStd] = deal(xHFlyStds(xHFlyMeans~=0), yHFlyStds(xHFlyMeans~=0));
+[xMin, yMin] = deal(xHFlyMins(xHFlyMeans~=0), yHFlyMins(xHFlyMeans~=0));
+[xMax, yMax] = deal(xHFlyMaxs(xHFlyMeans~=0), yHFlyMaxs(xHFlyMeans~=0));
+
+hFig = createPrintFig(7 * [1 1]);
+hScatter = errorbar(x, y, xMin - x, xMax - x, yMin - y, yMax - y, ...
+                    'LineStyle', 'none', 'CapSize', 0, 'LineWidth', 1, ...
+                    'Color', [1 1 1] * 0.4);
+% hScatter = errorbar(x, y, -xStd, xStd, -yStd, yStd, ...
+%                     'LineStyle', 'none', 'CapSize', 0, 'LineWidth', 1, ...
+%                     'Color', [1 1 1] * 0.2);
+hold on
+hScatter = scatter(x, y, 25, 'filled', 'MarkerFaceColor', [1 1 1] * 0.20);
+lsline
+if doBootstrap && numel(x) > 1
+    [bootstrapCorrCI, bootstrapCorrs] = bootci(10000, @corr, x(:), y(:));
+    legStr = sprintf('n = %d, corr = %.2f\n0.95-CI = [%.2f, %.2f]', ...
+                      numel(x), mean(bootstrapCorrs), bootstrapCorrCI);
+else
+    [corrCoef, p] = corr(x(:), y(:));
+    legStr = ['n = ' num2str(numel(x)) ...
+              ', \n corr = ' num2str(corrCoef, 2) ', \n p = ' num2str(p, 2)];
+end
+title(sprintf(legStr));
+xlabel('Tm9 fwhm vert. (deg)');
+ylabel('T5 fwhm vert. (deg)');
+% hLeg = legend(hScatter, legStr, 'Location', 'northoutside');
+% hLeg.Box = 'off';   
+prettifyAxes(gca);
+
+print(gcf, [figDir 'rfFwhmScatterFlyMeans-YOff-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');
+
+%% Check RF position distribution
+
+xH = vec(pos{1}(1, isAboveThreshold));
+yH = vec(pos{2}(1, isAboveThreshold));
+xV = vec(pos{1}(2, isAboveThreshold));
+yV = vec(pos{2}(2, isAboveThreshold));
+tmpInds = flyInds(isAboveThreshold);
+nValidFlies = numel(unique(tmpInds));
+iAx = 1;
+
+hFig = createPrintFig(7 * [ 1 2/3]);
+
+hHAx(1) = subplot(1, 2, 1);
+hHist(1) = histogram(xH - yH);
+hHAx(2) = subplot(1, 2, 2);
+hHist(2) = histogram(xV - yV);
+
+[hHist.EdgeColor] = deal('w');
+[hHist.FaceColor] = deal([1 1 1] * 0.0);
+
+title(hHAx(1), {'Hor. RF Pos.', '(Tm9 - T5)'});
+title(hHAx(2), {'Hor. RF Pos.', '(Tm9 - T5)'});
+
+xlabel(hHAx(1), 'RF Pos. Diff. (deg)');
+
+arrayfun(@prettifyAxes, hHAx);
+setFontForThesis(hHAx, hFig);
+
+print(gcf, [figDir 'rfPosDiffScatter-Off-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');
+
+%% Extract the parameters for the analysis of spatial RFs.
+% clear amp pos width rsquare responseIndex
+tuningMethod = 'median';
+StimTableTmp = StimTableNeuronStim{2};
+tuningInd = strcmp(tuningMethod, StimTableTmp.nonParamTuning{1}(1).tuningMethods);
+jBar = 1;
+for iStim = 2:5
+    dummyTable = StimTableNeuronStim{1, iStim};
+    for jChannel = 1: 2
+        responseIndex{jChannel}(jBar, :) = cell2mat(cellfun(@(x, y) x(jChannel, 1: end - 1*0), dummyTable.responseIndex, 'uni', 0)');        
+        % Hard coded getting rid of background ROI, but this should also be
+        % imporved to handle the roiMeta field to discard invalid ROIs
+        % (rois with NaNs)
+        tuningCurves{jChannel}(jBar, :) = cellfun(@(x) squeeze(x(jChannel).tc{tuningInd}(:, 1: end - 1, :)), ...
+                                                  dummyTable.nonParamTuning, 'uni', 0);
+    end
+    jBar = jBar + 1;
+end
+%% Plot tuning curves.
+for iStim = 1: 4
+    tm9AllTCs{iStim} = cat(1, tuningCurves{1}{iStim, :});
+    t5AllTCs{iStim} = cat(1, tuningCurves{2}{iStim, :});
+end
+% Align to fit
+% pos is already original parameter b1 as 2*(b1 - 1)
+
+lagsCell = cellfun(@(x, y) round(x/2+1 - size(y, 2) / 2), ...
+                   mat2cell(pos{1}, [1 1 1 1])', tm9AllTCs, 'uni', 0);
+Tm9BarRFs = cellfun(@(x, y) alignRFsToAbsMax(x, 0, y'), tm9AllTCs, lagsCell, 'uni', 0);
+
+lagsCell = cellfun(@(x, y) round(x/2+1 - size(y, 2) / 2), ...
+                   mat2cell(pos{2}, [1 1 1 1])', t5AllTCs, 'uni', 0);
+T5BarRFs = cellfun(@(x, y) alignRFsToAbsMax(x, 0, y'), t5AllTCs, lagsCell, 'uni', 0);
+
+
+%% Blank epoch was already deleted, now delete vertical bars on edges of screen.
+Tm9BarRFs{2}(:, end-1: end) = [];
+Tm9BarRFs{4}(:, end-1: end) = [];
+Tm9BarRFs{2}(:, 1: 2) = [];
+Tm9BarRFs{4}(:, 1: 2) = [];
+
+T5BarRFs{2}(:, end-1: end) = [];
+T5BarRFs{4}(:, end-1: end) = [];
+T5BarRFs{2}(:, 1: 2) = [];
+T5BarRFs{4}(:, 1: 2) = [];
+%% Plot only analyzed OFF RFs.
+minRespInd = min(responseIndexOff);
+thresholdRSquare = 0.25;
+minRSquare = min(minOffRSquare);
+isAboveThreshold = minRSquare > thresholdRSquare & minRespInd > 0.3;
+Tm9BarRFs = cellfun(@(x) x(isAboveThreshold, :), Tm9BarRFs, 'uni', 0);
+T5BarRFs = cellfun(@(x) x(isAboveThreshold, :), T5BarRFs, 'uni', 0);
+%%
+% Plot lasagna plots, separate figures for off and on, and reapeated panel
+% is to have a colorbar that does not distort existing axis of interest.
+hFig = createPrintFig(15 * [ 1 1/4]);
+hSubAx = plotLasagnaTracesFromCell(gcf, [Tm9BarRFs, Tm9BarRFs(1)], brewermap(6, 'Paired'));
+colorbar
+arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx([1, 3], 1));
+arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'off'), hSubAx(2: end, 2));
+arrayfun(@prettifyAxes, hSubAx)
+arrayfun(@offsetAxes, hSubAx)
+arrayfun(@(x) set(x, 'Color', 'none'), hSubAx);
+
+% set(hSubAx, 'Visible', 'off')
+linkaxes(hSubAx(1: 2, 1), 'y')
+linkaxes(hSubAx(3: 4, 1), 'y')
+hFig.Color = 'none';
+figFileName = ['Tm9-RFcurves-OFF-XY-lasagna'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+%%
+hFig = createPrintFig(15 * [ 1 1/4]);
+hSubAx = plotLasagnaTracesFromCell(gcf, [T5BarRFs, T5BarRFs(1)], brewermap(6, 'Paired'));
+colorbar
+arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx([1, 3], 1));
+arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'off'), hSubAx(2: end, 2));
+arrayfun(@prettifyAxes, hSubAx)
+arrayfun(@offsetAxes, hSubAx)
+arrayfun(@(x) set(x, 'Color', 'none'), hSubAx);
+
+% set(hSubAx, 'Visible', 'off')
+linkaxes(hSubAx(1: 2, 1), 'y')
+linkaxes(hSubAx(3: 4, 1), 'y')
+hFig.Color = 'none';
+figFileName = ['T5-RFcurves-OFF-XY-lasagna'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+
+%% Export data for figure 6b-e
+
+receptiveFields = [Tm9BarRFs; T5BarRFs];
+dataFileName = [figDir filesep 'Fig6.xls'];
+cellTypeStr = {'Tm9', 'T5'};
+sheetNamesCell = {'bc', 'de'};
+
+stimNameStrCell = {'_OFF_horizontal_position_deg', '_OFF_vertical_position_deg', ...
+                   '_ON_horizontal_position_deg', '_ON_vertical_position_deg'};
+               
+for iCellType = 1: numel(cellTypeStr)
+    startRow = 1;
+    sheetName = ['Fig6' sheetNamesCell{iCellType}];
+    for iStim = 1: 4
+        parameterLabel = [cellTypeStr{iCellType} stimNameStrCell{iStim}];
+        DataTable = createRFsTable(receptiveFields{iCellType, iStim}, parameterLabel);
+        writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
+        startRow = startRow + size(DataTable, 1) + 2;
+    end
+end
+%% Tests if T5 ON response is 0 for reviewer #2
+
+xH = vec(amp{1}(3, isAboveThreshold));
+yH = vec(amp{2}(3, isAboveThreshold));
+xV = vec(amp{1}(4, isAboveThreshold));
+yV = vec(amp{2}(4, isAboveThreshold));
+% use minimum for negative ON response.
+minPeaks = cellfun(@(x) min(x, [], 2), T5BarRFs(3: 4), 'uni', 0);
+tcStd = cellfun(@(x) std(x, [], 2), T5BarRFs(3: 4), 'uni', 0);
+peakStdRatio = cellfun(@(x, y) x ./ y, minPeaks, tcStd, 'uni', 0);
+
+% figure
+% for iAx = 1: 2
+%     subplot(1, 2, iAx);
+%     histogram(peakStdRatio{iAx})
+% end
+ % check if the fitted amplitude is higher than the tuning curve values
+ % (noise), signrank uses x - y, so tc - amp shuould be positive for a
+ % smaller amplitude (larger negative peak), so we need right tail test.
+[pValCell, isSignificant] = cellfun(@(x, y) arrayfun(@(z) signrank(x(z, :), y(z), 'tail', 'right'), 1: numel(y)), ...
+                   T5BarRFs(3: 4), {yH, yV}, 'uni', 0);
+
+proportionSignificant = cellfun(@(x) sum(x) / numel(x), isSignificant);
+
+clear hAx
+hFig = createPrintFig(15 * [3/4 1]);
+hAx(1) = subplot(3, 2, 1);
+histogram(yH);
+hAx(2) = subplot(3, 2, 2);
+histogram(yV);
+
+
+for iAx = 1: 2
+    hAx(2 + iAx) = subplot(3, 2, 2  * iAx + 1);
+    histogram(pValCell{iAx}, 'NumBins', 100);
+end
+hAx(5) = subplot(3, 2, [4, 6]);
+hBar = bar([proportionSignificant; 1 - proportionSignificant]', 'stacked');
+hLeg = legend(hAx(5), {'p < 0.05', 'p >= 0.05'}, 'Location', 'best');
+hAx(5).XTickLabels = {'Hor.',  'Vert.'};
+ylim([0 1])
+
+suptitle('T5 Responses to ON bars')
+
+arrayfun(@prettifyAxes, hAx);
+setFontForThesis(hAx, hFig);
+
+xLabelCell = {'Amplitude Hor. (\Delta F / F_0)', ...
+              'Amplitude Vert. (\Delta F / F_0)', ...
+              '', 'p-Value', 'Bar orientation'};
+yLabelCell = {'No. ROIs', '', 'No. ROIs', 'No. ROIs', 'Proportion of ROIs'};
+
+arrayfun(@(x, y) xlabel(x, y), hAx, xLabelCell);
+arrayfun(@(x, y) ylabel(x, y), hAx, yLabelCell);
+
+figFileName = ['T5-ON-Amp-XY-stats'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+
+%%
+
+function tcRoiCell = getTuningCurvePerROI(tc, tuningMethods, selectedTuningMethod)
+%% Expose single Roi tuning maps.
+% The epochs are varying first in the spacing then in the temporal
+% frequency, so tc consists of 1 x nTFreqs and each element is of shape
+% 1 x nRois x nSFreqs.
+% Here we first extract per Roi the nSFreqs tuning, as column vector, and
+% concatenate the nTFreqs column vectors arising this way to get a tuning
+% matrix of nSFreqs x nTFreqs.
+% tc1F = tc(:, strcmp(tuningMethods, '1F_amplitude'));
+tc = tc(:, strcmp(tuningMethods, selectedTuningMethod));
+tcRoiCell = cell(1, size(tc{1}, 2));
+
+for iRoi = 1: size(tc{1}, 2)
+    tcRoiCell{iRoi} = cell2mat(cellfun(@(x) squeeze(x(:, iRoi, :)), tc, 'uni', 0)');
+end
+
+if  0
+% Try to reshape the tc into one tc array for each roi. is equivalent but
+% maybe less clear. First concatenate along stim dimension to get an array
+% of dims 1 x nRois x (nSFreqs x nTFreqs), squeeze it to get an array of
+% dims nRois x (nSFreqs x nTFreqs), and reshape it into an array of dims
+% nRois x nSFreqs x nTFreqs. Now we slice it along the first dimension into
+% a cell array of nRois x 1, each cell being nSFreqs x nTFreqs.
+tcRoiCell = mat2cell(reshape(squeeze(cat(3, tc{:})), sizeTC(2), sizeTC(3), []), ...
+                    ones(1, sizeTC(2)), stimSizeTC(1), stimSizeTC(2));
+end                
+
+end

+ 350 - 0
figures/MakeRFNoiseComparison_Fig2ef.m

@@ -0,0 +1,350 @@
+%% Make fig 2 e, f
+close all;
+clear all;
+clc;
+
+%% Define search path.
+parentDir = fileparts(fileparts(mfilename('fullpath')));
+saveFolder = [parentDir filesep 'data'];
+load([saveFolder filesep 'Fig2-Tm9-Bars-Noise-processedTableWithBackgroundOnOffNoiseCropped'], 'ProcessedTable');
+
+%%
+figDir = [parentDir filesep 'figures' filesep 'Fig2' filesep];
+if ~exist(figDir, 'dir')
+    mkdir(figDir); 
+end
+
+%% Define search path.
+datasetInd = 3;
+[~, ~, stimSetCellStr] = parseDatasetVars(datasetInd);
+
+%%
+ProcessedTable = ProcessedTable(ProcessedTable.typeOfProblem == 0, :);
+% Check only bars.
+nStims = 8; % 8 for fff, bars, and noise; 5 for fff and bars.
+stimSetCellStr = stimSetCellStr(1: nStims);
+%% Split table into cell types.
+tokens = arrayfun(@(x)  regexp(x, '.*Tm(\d).*', 'tokens'), ProcessedTable.timeSeriesPath);
+cellTypeNumArray = cellfun(@(x) str2num(x{1}{1}), tokens);
+cellTypeNum = unique(cellTypeNumArray);
+for iNeuron = 1: numel(cellTypeNum)
+    StimTableNeuronCell{iNeuron} = ProcessedTable(cellTypeNumArray == cellTypeNum(iNeuron), :);
+end
+%%
+stimInd = 1;
+nStims = numel(stimSetCellStr);
+for iNeuron = 1: numel(cellTypeNum)
+    for jStim = 1: nStims
+        StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
+    end
+end
+
+for iFly = 1: numel(unique(ProcessedTable.flyInd))
+    flyInds = ProcessedTable.flyInd == iFly;
+    flyStims = ProcessedTable(flyInds, :).stimParamFileName;
+    % Find if stim is within the chosen subset.
+    flyStimInds{iFly} = (cellfun(@(x) find(strcmp(x, stimSetCellStr)), flyStims, 'uni', 0));
+    % If it contains stimuli, check that it includes the 5 chosen.
+    if ~isempty(flyStimInds{iFly})
+        setDiff{iFly} = setdiff(1: nStims, [flyStimInds{iFly}{:}]);
+    else
+        setDiff{iFly} = -1;
+    end
+    % For flies with all chosen stimuli order the table accordingly.
+    if isempty(setDiff{iFly})
+        sortFlyIndsTemp = find(flyInds);
+        [~, sortInd] = sort([flyStimInds{iFly}{:}]);
+        sortFlyInds{iFly} = sortFlyIndsTemp(sortInd);
+    end
+end
+
+StimOrderedTable = ProcessedTable(cat(1, sortFlyInds{:}), :);
+%% Split table into cell types.
+
+tokens = arrayfun(@(x)  regexp(x, '.*Tm(\d).*', 'tokens'), StimOrderedTable.timeSeriesPath);
+cellTypeNumArray = cellfun(@(x) str2num(x{1}{1}), tokens);
+cellTypeNum = unique(cellTypeNumArray);
+for iNeuron = 1: numel(cellTypeNum)
+    StimTableNeuronCell{iNeuron} = StimOrderedTable(cellTypeNumArray == cellTypeNum(iNeuron), :);
+end
+%%
+nStims = numel(stimSetCellStr);
+for iNeuron = 1: numel(cellTypeNum)
+    for jStim = 1: nStims
+        StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
+    end
+end
+
+OriginalStimTableNeuronStim = StimTableNeuronStim;
+
+%% Delete invalid rois for all flies and stimuli.
+
+bgRois = cellfun(@(y) cellfun(@(x) x.backgroundRois, y.roiMeta), StimTableNeuronStim, 'uni', 0);
+%Test if all background rois are the same for all stimuli.
+assert(isequal(bgRois{:}));
+invalidRois = zeros(size(StimTableNeuronStim{1}, 1), size(StimTableNeuronStim, 2));
+for iFly = 1: size(StimTableNeuronStim{1}, 1)
+    for jStim = 1: size(StimTableNeuronStim, 2)
+        invalidInds = StimTableNeuronStim{jStim}.roiMeta{iFly}.invalidRois;
+        if ~isempty(invalidInds)
+            invalidRois(iFly, jStim) = invalidInds;
+        end
+    end
+end
+
+roisToDeletePerFly = sum(unique(invalidRois', 'rows'))';
+
+%% Extract the parameters for the analysis of spatial RFs.
+jBar = 1;
+for iStim = 2:5
+    dummyTable = StimTableNeuronStim{iStim};
+    paramsCell = cellfun(@(x) x(1: end - 1*0), dummyTable.paramTuning, 'uni', 0);
+    responseIndex = cellfun(@(x, y) x(1: end - 1*0), dummyTable.responseIndex, 'uni', 0);
+    for jFly = find(roisToDeletePerFly ~= 0)
+        if invalidRois(jFly, iStim) == 0        
+            paramsCell{jFly}(roisToDeletePerFly(jFly) - 1*0) = [];
+            responseIndex{jFly}(roisToDeletePerFly(jFly) - 1*0) = [];
+        end
+    end
+    amp(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.a1, x, 'uni', 1), paramsCell, 'uni', 0)');
+    pos(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.b1, x, 'uni', 1), paramsCell, 'uni', 0)');
+    width(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.c1, x, 'uni', 1), paramsCell, 'uni', 0)');
+    rsquare(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.gof.rsquare, x, 'uni', 1), paramsCell, 'uni', 0)');
+    jBar = jBar + 1;
+end
+
+
+for iStim = [7, 8]
+    dummyTable = StimTableNeuronStim{iStim};
+    paramsCell = cellfun(@(x) x(1: end - 1*0), dummyTable.paramTuning, 'uni', 0);
+%     responseIndex = cellfun(@(x, y) x(1: end - 1), dummyTable.responseIndex, 'uni', 0);                    
+    for jFly = find(roisToDeletePerFly ~= 0)
+        if invalidRois(jFly, iStim) == 0
+            paramsCell{jFly}(roisToDeletePerFly(jFly) - 1*0) = [];
+%             responseIndex{jFly}(roisToDeletePerFly(jFly) - 1*0) = [];
+        end
+    end
+    amp(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.a1, x, 'uni', 1), paramsCell, 'uni', 0)');
+    pos(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.b1, x, 'uni', 1), paramsCell, 'uni', 0)');
+    width(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.c1, x, 'uni', 1), paramsCell, 'uni', 0)');
+    rsquare(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.gof.rsquare, x, 'uni', 1), paramsCell, 'uni', 0)');
+    jBar = jBar + 1;
+end
+
+pos(1: 4, :) = 2 * (pos(1: 4, :) - 1);
+pos(5: 6, :) = 5 * pos(5: 6, :);
+width = (2 * sqrt(log(2))) * width;
+width(1: 4, :) = 2 * width(1: 4, :);
+width(5: 6, :) = 5 * width(5: 6, :);
+validInds = all(rsquare([1, 2, 5, 6], :) > 0.5) & all(rsquare([3, 4], :) > 0.1);
+%%
+flyInds = repelem(1:size(StimTableNeuronStim{1}),cellfun(@numel, StimTableNeuronStim{1}.responseIndex));
+flyInds(roisToDeletePerFly(1)) = [];
+flyInds(validInds)
+numel(unique(flyInds(validInds)))
+%%
+clear allRFs
+for iStim = 2: 5
+    tuningMethods = StimTableNeuronStim{iStim}.nonParamTuning{1}.tuningMethods;
+    tuningMethodInd = strcmp(tuningMethods, 'extreme');
+    tcCellArray = cellfun(@(x) squeeze(x.tc{tuningMethodInd}(:, 1: end - 1, :)), StimTableNeuronStim{iStim}.nonParamTuning, 'uni', 0);
+    allRFs{iStim - 1} = cell2mat(tcCellArray);
+end
+%%
+noiseRFs{1} = getSpaceFilterArrayFromTable(StimTableNeuronStim{7}); % Az is horizontal
+noiseRFs{2} = getSpaceFilterArrayFromTable(StimTableNeuronStim{8}); % El is vertical
+%%
+% Delete invalid roi that shifts all images by one;                                                       
+barRFs = cellfun(@(x) x(setdiff(1: length(allRFs{1}), roisToDeletePerFly(1)), 1: end), ...
+                 allRFs, 'uni', 0);
+noiseRFs{1}(roisToDeletePerFly(1), :) = [];
+
+%%             
+% The screen is about 48-60 deg, we center the RFs at about the center 24
+% deg, using the fit positions for higher presicion.
+lagsBars = round((pos(1: 4, :) - 25) / 2);
+for iBar = 1: numel(barRFs)
+    barRFs{iBar} = alignRFsToAbsMax(barRFs{iBar}, 0, lagsBars(iBar, :)');
+end
+
+lagsNoise = round((pos(5: 6, :) - 30) / 5);
+for iBar = 1: numel(noiseRFs)
+    noiseRFs{iBar} = alignRFsToAbsMax(noiseRFs{iBar}, 0, lagsNoise(iBar, :));
+end
+
+%% Blank epoch was already deleted, now delete vertical bars on edges of screen.
+barRFs{2}(:, end-1: end) = [];
+barRFs{4}(:, end-1: end) = [];
+barRFs{2}(:, 1: 2) = [];
+barRFs{4}(:, 1: 2) = [];
+%% Fig2e
+hFig = figure;
+splitSign = 1;
+doSmoothing = 1;
+indsToPlot = validInds;
+subplot_tight(1, 3, 1);
+plotRFEnsemble(noiseRFs{1}, noiseRFs{2}, splitSign, doSmoothing, indsToPlot);
+title('White noise spatial filters')
+% Plot bar RFs. OFF.
+subplot_tight(1, 3, 2);
+plotRFEnsemble(barRFs{1}, barRFs{2}, splitSign, doSmoothing, indsToPlot);
+title('OFF flashing bar filters')
+% Plot bar RFs. ON.
+subplot_tight(1, 3, 3);
+plotRFEnsemble(barRFs{3}, barRFs{4}, splitSign, doSmoothing, indsToPlot);
+title('ON flashing bar filters')
+suptitle('Cartesian product (X \times Y) of one dimensional tuning curves');
+
+figFileName = ['RFIm-ON-OFF-Noise'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+% print(hFig, [figDir figFileName '.png'], '-dpng', '-r
+%% Write data for fig2e
+dataFileName = [figDir filesep 'Fig2.xls'];
+DataTable = createRFsTable(noiseRFs{1}(indsToPlot, :), 'Noise_horizontal_position_deg');
+% Correct spacing from 2 deg for bars to 5 deg for noise.
+DataTable(:, 1).Noise_horizontal_position_deg = 5 / 2 * DataTable(:, 1).Noise_horizontal_position_deg;
+writetable(DataTable, dataFileName, 'Sheet', 'Fig2e')
+DataTable = createRFsTable(noiseRFs{1}(indsToPlot, :), 'Noise_vertical_position_deg');
+% Correct spacing from 2 deg for bars to 5 deg for noise.
+DataTable(:, 1).Noise_vertical_position_deg = 5 / 2 * DataTable(:, 1).Noise_vertical_position_deg;
+writetable(DataTable, dataFileName, 'Sheet', 'Fig2e', 'Range', 'A15')
+
+DataTable = createRFsTable(barRFs{1}(indsToPlot, :), 'OFF_horizontal_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig2e', 'Range', 'A29')
+DataTable = createRFsTable(barRFs{2}(indsToPlot, :), 'OFF_vertical_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig2e', 'Range', 'A56')
+DataTable = createRFsTable(barRFs{3}(indsToPlot, :), 'ON_horizontal_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig2e', 'Range', 'A83')
+DataTable = createRFsTable(barRFs{4}(indsToPlot, :), 'ON_vertical_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig2e', 'Range', 'A110')
+
+%%
+offX = barRFs{1}(validInds, :);
+offY = barRFs{2}(validInds, :);
+onX = barRFs{3}(validInds, :);
+onY = barRFs{4}(validInds, :);
+lnX = noiseRFs{1}(validInds, :);
+lnY = noiseRFs{2}(validInds, :);
+
+lnX = interp1(linspace(2.5,60, 12), lnX', linspace(2.5,60, 25))';
+lnY = interp1(linspace(2.5,60, 12), lnY', linspace(2.5,60, 25))';
+% Realign noise to center now at 2 deg precision.
+[~, peakInd] = max(abs(lnX), [], 2);
+lnX = alignRFsToAbsMax(lnX, 0, peakInd - round(length(lnX) / 2) - 1);
+[~, peakInd] = max(abs(lnY), [], 2);
+lnY = alignRFsToAbsMax(lnY, 0, peakInd - round(length(lnY) / 2) - 1);
+
+maxNorm = @(x, dim) x ./ max(abs(x), [], dim);
+getRow = @(x, nRow) x(nRow, :);
+normRFs = cellfun(@(x) cell2mat(arrayfun(@(y) smooth(getRow(maxNorm(x,2), y))', 1: size(x, 1), 'uni', 0)'), ...
+                  {offX, offY, onX, onY, lnX, lnY}, 'uni', 0);
+normRFs = cellfun(@(x) maxNorm(x,2), normRFs, 'uni', 0);
+
+hFig = createPrintFig(15 * [ 1 1 / 3]);
+for iRFType = 1:3
+    subplot_tight(1, 3, iRFType);
+    imshow(getRFIm(mean(normRFs{2*iRFType}), mean(normRFs{2*iRFType - 1})))
+end
+
+figFileName = ['RFIm-ON-OFF-Noise-mean'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+
+hFig = createPrintFig(15 * [ 1 2/3]);
+hSubAx = plotLasagnaTracesFromCell(gcf, normRFs, brewermap(6, 'Paired'));
+ylim(hSubAx(:, 1), [-1 1]);
+arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx(1, :))
+cellfun(@(x) caxis(ZeroCenteredBounds(x, [0 100])), normRFs, 'uni', 0)
+arrayfun(@prettifyAxes, hSubAx)
+arrayfun(@offsetAxes, hSubAx)
+set(hSubAx, 'Visible', 'off')
+figFileName = ['RFcurves-ON-OFF-Noise-XY-lasagna'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+
+%%
+hFig = createPrintFig(15 * [1 1/2]);
+pairsToAverage = {[1, 2], [3, 4], [5, 6]};
+meanWidths = cellfun(@(x) mean(width(x, validInds)), pairsToAverage, 'uni', 0);
+meanWidths = cellfun(@(x) mean(width(x, validInds)), pairsToAverage, 'uni', 0);
+
+beanPlot(meanWidths, {'off', 'on', 'noise'}, flipud(brewermap(3, 'Set1')), 1, gca, [0-eps 50], 'vertical', false, true)
+figFileName = ['RF-ON-OFF-Noise-FWHM-Bean'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+%% Save data for Fig 2f
+xlswrite(dataFileName, {'FHWM-OFF-bars'}, 'Fig2f', 'A1')
+xlswrite(dataFileName, meanWidths{1}, 'Fig2f', 'B1')
+xlswrite(dataFileName, {'FHWM-ON-bars'}, 'Fig2f', 'A2')
+xlswrite(dataFileName, meanWidths{2}, 'Fig2f', 'B2')
+xlswrite(dataFileName, {'FHWM-Noise-bars'}, 'Fig2f', 'A3')
+xlswrite(dataFileName, meanWidths{3}, 'Fig2f', 'B3')
+%% Permutation test for mean RF width difference.
+[pValOFFvsON, obsDiff] = permutationTest(meanWidths{1},meanWidths{2}, 100000)
+[pValONvsNoise, obsDiff] = permutationTest(meanWidths{2},meanWidths{3}, 100000)
+[pValOFFvsNoise, obsDiff] = permutationTest(meanWidths{1},meanWidths{3}, 100000)
+%%
+hFig = createPrintFig(15 * [1/2 1/2]);
+
+meanAmps = cellfun(@(x) mean(amp(x, validInds)), pairsToAverage(1:2), 'uni', 0);
+
+beanPlot(meanAmps, {'off', 'on'}, flipud(brewermap(3, 'Set1')), 1, gca,'unbounded', 'vertical', false, true)
+figFileName = ['RF-ON-OFF-Noise-Amp-Bean'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+
+%% Save data for Fig 2f
+xlswrite(dataFileName, {'Amp-OFF-bars'}, 'Fig2f', 'A5')
+xlswrite(dataFileName, meanAmps{1}, 'Fig2f', 'B5')
+xlswrite(dataFileName, {'Amp-ON-bars'}, 'Fig2f', 'A6')
+xlswrite(dataFileName, meanAmps{2}, 'Fig2f', 'B6')
+%%
+% pairsToPlot = {[1 5], [2 6], [3 5], [4 6]};
+pairsToPlot = {[1 2], [3 4], [1 4], [1 3], [2 4], [2 3]};
+labelsCell = {'X OFF', 'Y OFF', 'X ON', 'Y ON', 'X Noise', 'Y Noise'};
+indsToPlot = validInds;
+% fwhmArray = [mean(width(1: 2, :)); mean(width(3: 4, :)); mean(width(5: 6, :))]';
+fwhmArray = width';
+rSquare = rsquare';
+% hFig = createFullScreenFig;
+% set(hFig, 'Units', 'Centimeters')
+% set(hFig, 'Position', [hFig.Position(1:2) + 10 14 14]);
+hFig = createPrintFig(15 * [ 1 2/3]);
+
+
+for iPair = 1: numel(pairsToPlot)
+    hScattAx(iPair) = subplot(2, 3, iPair); hold on;
+    xInd = pairsToPlot{iPair}(1);
+    yInd = pairsToPlot{iPair}(2);
+    scatter(fwhmArray(indsToPlot, xInd), fwhmArray(indsToPlot, yInd), [], ...
+            min(rSquare(indsToPlot, xInd), rSquare(indsToPlot, yInd)), 'filled');  
+    [correlation, pValue] = corr(fwhmArray(indsToPlot, xInd), fwhmArray(indsToPlot, yInd), 'type', 'Pearson' );
+    x = fwhmArray(indsToPlot, xInd);
+    [b, bInt]= regress(fwhmArray(indsToPlot, yInd), [ones(length(x), 1) x]);
+    hFit = plot(x, b(1) + x*b(2), 'Color', [1 1 1] * 0.5);
+    maxVal = max([fwhmArray(indsToPlot, xInd); fwhmArray(indsToPlot, yInd)]);
+    maxVal = quantile(vec(fwhmArray(indsToPlot, :)), 1);
+    plot([0 maxVal], [0 maxVal], 'Color', [1 1 1] * 0.5);
+    prettifyAxes(gca);
+    set(gca, 'XLim', [0 maxVal], 'YLim', [0 maxVal], 'CLim', [0 1]);
+    axis square;
+    xlabel(labelsCell{xInd});
+    ylabel(labelsCell{yInd});
+    title({['\rho = ' num2str(correlation, 2) ', p = ' num2str(pValue, 2) ],...
+           ['m = ' num2str(b(2), 2) ' (' num2str(bInt(2, 1), 2) ', ' num2str(bInt(2, 2), 2) ')']})
+    setFavoriteColormap
+end
+setFontForThesis(hScattAx, gcf);
+
+arrayfun(@prettifyAxes, hScattAx)
+arrayfun(@offsetAxes, hScattAx)
+% set(hSubAx, 'Visible', 'off')
+figFileName = ['BarRFs-Corr-Pearson'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+%%
+offWidth = mean(width(1: 2, indsToPlot), 1);
+onWidth = mean(width(3: 4, indsToPlot), 1);
+[corrWidth, pCorrWidth] = corr(offWidth(:), onWidth(:))
+
+%%
+offWidth = width(1: 2, indsToPlot);
+onWidth = width(3: 4, indsToPlot);
+[corrWidth, pCorrWidth] = corr(offWidth(:), onWidth(:))
+%%
+[corrWidth, pCorrWidth] = corr(width(1:4, indsToPlot)', 'type','Spearman' )

+ 334 - 0
figures/MakeTm4Tm9RFs_fig3.m

@@ -0,0 +1,334 @@
+ccc
+%% Define search path.
+parentDir = fileparts(fileparts(mfilename('fullpath')));
+saveFolder = [parentDir filesep 'data'];
+load([saveFolder filesep 'Fig3-data_Tm9GCaMP6f-Tm4jRGECO1a_ONOFFPaperNew']);
+
+%%
+figDir = [parentDir filesep 'figures' filesep 'Fig3' filesep];
+if ~exist(figDir, 'dir'); mkdir(figDir); end
+
+%% Genereta figure 3b
+onOffPerCell{1} = cell2mat(cellfun(@(x) cell2mat(cellfun(@(y) y, x(:, 1), 'uni', 0)), onoffData, 'Uni', 0)')';
+onOffPerCell{2} = cell2mat(cellfun(@(x) cell2mat(cellfun(@(y) y, x(:, 2), 'uni', 0)), onoffData, 'Uni', 0)')';
+timePoints = (0: size(onOffPerCell{1}, 2)-1) /10;
+
+plotONOFFTracesFromCell(onOffPerCell, {'Tm4', 'Tm9'}, [figDir filesep], brewermap(2, 'PrGn'));
+
+%% Save data from plot.
+dataFileName = [figDir filesep 'Fig3.xls'];
+DataTable = createRFsTable(onOffPerCell{1}, 'Time_seconds_Tm4');
+DataTable(:, 1).Time_seconds_Tm4 = timePoints(:);
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3b', 'Range', 'A1');
+
+DataTable = createRFsTable(onOffPerCell{2}, 'Time_seconds_Tm9');
+DataTable(:, 1).Time_seconds_Tm9 = timePoints(:);
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3b', 'Range', ['A' num2str(size(onOffPerCell{1}, 2) + 3)])
+
+%% Reanalyze data of dual imaging.
+ccc
+
+parentDir = fileparts(fileparts(mfilename('fullpath')));
+saveFolder = [parentDir filesep 'data'];
+load([saveFolder filesep 'Fig3-data_Tm9GCaMP6f-Tm4jRGECO1a_All']);
+
+%%
+figDir = [parentDir filesep 'figures' filesep 'Fig3' filesep];
+if ~exist(figDir, 'dir'); mkdir(figDir); end
+%%
+[nFlies, nLayers, nStims] = size(dataEdges);
+% New functions use cell arrays.
+dataEdges = mat2cell(dataEdges, ones(1, nFlies), ones(1, nLayers), ones(1, nStims));
+% Get tuning curves.
+tuningStat = 'maxAbs';
+% tuningStat = 'median';
+barRFs = collectRFs(dataEdges, tuningStat);
+
+% Remove blanck epoch.
+barRFs = cellfun(@(x) x(:, 2: end), barRFs, 'Uni', 0);
+barRFs([2, 4], :) = cellfun(@(x) x(:, 3: end - 2), barRFs([2, 4],  :), 'Uni', 0);
+
+%%
+
+fitStructOff = cellfun(@(x) fitGaussian1(1: size(x, 2), x, 1), barRFs(1:2, :), 'Uni', 0);
+fitStructOn = cellfun(@(x) fitGaussian1(1: size(x, 2), x, -1), barRFs(3:4, :), 'Uni', 0);
+fitStructBars = [fitStructOff; fitStructOn];
+
+%% Extract the parameters for the analysis of spatial RFs.
+for iChannel = 1: 2
+    amp{iChannel} = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.a1, x, 'uni', 1), ...
+                                        fitStructBars(:, iChannel), 'uni', 0));
+    pos{iChannel} = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.b1, x, 'uni', 1), ...
+                                        fitStructBars(:, iChannel), 'uni', 0));
+    width{iChannel} = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.c1, x, 'uni', 1), ...
+                                          fitStructBars(:, iChannel), 'uni', 0));
+    rsquare{iChannel} = cell2mat(cellfun(@(x) arrayfun(@(y) y.gof.rsquare, x, 'uni', 1), ...
+                                            fitStructBars(:, iChannel), 'uni', 0));
+end
+pos = cellfun(@(x) 2 * (x - 1), pos, 'uni', 0);
+width = cellfun(@(x) 2 * (2 * sqrt(log(2))) * x, width, 'uni', 0);
+
+%%
+% barRFs = cellfun(@(x) alignRFsToAbsMax(x, 0), barRFs, 'uni', 0);
+
+% The screen is about 48-60 deg, we center the RFs at about the center 24
+% deg, using the fit positions for higher presicion.
+for iChannel = 1: 2
+    lagsBars = round((pos{iChannel}(1: 4, :) - 25) / 2);
+    for iBar = 1: size(barRFs, 1)
+        barRFs{iBar, iChannel} = alignRFsToAbsMax(barRFs{iBar, iChannel}, 0, lagsBars(iBar, :)');
+    end
+end
+%% Use all data for fig3c-d
+splitSign = 1;
+doSmoothing = 1;
+rSquareCutoff = -20;
+plotImage = 1;
+hue = [];
+useSubplots = 1;
+
+for iChannel = 1: 2
+    rSquareArray = rsquare{iChannel}';
+    % Check good fit for all stim.
+    [rfInds{1: 2, iChannel}] =  deal(all(rSquareArray(:, [1: 2]) >= rSquareCutoff, 2));
+    [rfInds{3: 4, iChannel}] =  deal(all([rSquareArray(:, [1:2]) >= rSquareCutoff ... 
+                                          rSquareArray(:, [3:4]) >= -20], 2));
+end
+%% Figure 3c, examples extracted from here.
+hFig = createPrintFig(10 * [ 1 1]);
+for iChannel = 1: 2
+    subplot(2, 2, 2 * iChannel - 1);
+    plotRFEnsemble(barRFs{2, iChannel}, barRFs{1, iChannel}, ...
+                   splitSign, doSmoothing, rfInds{1, iChannel}, plotImage, hue, useSubplots)
+    subplot(2, 2, 2 * iChannel);
+    plotRFEnsemble(barRFs{4, iChannel}, barRFs{3, iChannel}, ...
+                   splitSign, doSmoothing, rfInds{3, iChannel}, plotImage, hue, useSubplots)
+end
+figFileName = ['RFIm-ON-OFF-Noise-raw-ON-norm'];    
+print(hFig, [figDir filesep figFileName '.pdf'], '-dpdf')
+%% For the mean images in fig. 3d
+filteredRFs = cellfun(@(x, y) x(y, :), barRFs, rfInds, 'uni', 0);
+
+maxNorm = @(x, dim) x ./ max(abs(x), [], dim);
+getRow = @(x, nRow) x(nRow, :);
+% normRFs = cellfun(@(x) cell2mat(arrayfun(@(y) smooth(getRow(maxNorm(x,2), y))', 1: size(x, 1), 'uni', 0)'), ...
+%                   barRFs, 'uni', 0);
+normRFs = cellfun(@(x) cell2mat(arrayfun(@(y) smooth(getRow(x, y))', 1: size(x, 1), 'uni', 0)'), ...
+                  filteredRFs, 'uni', 0);
+normRFs = cellfun(@(x) maxNorm(x,2), normRFs, 'uni', 0);
+
+% Maximum of the mean tuning curve to normalize image accross stimuli.
+maxPerStim = max(cellfun(@(x) max(abs(mean(x))), filteredRFs, 'uni', 1), [], 2);
+maxPerPolarity = max(reshape(maxPerStim, [2 2]));
+hFig = createPrintFig(15 * [ 1 1 / 3]);
+for iRFType = 1:4
+    subplot_tight(1, 4, iRFType);
+%     imshow(getRFIm(mean(normRFs{2*iRFType}), mean(normRFs{2*iRFType - 1})))
+    imshow(getRFIm(mean(filteredRFs{2*iRFType}), mean(filteredRFs{2*iRFType - 1}), ...
+           maxPerPolarity(2 - mod(iRFType, 2))))
+end
+%%
+figFileName = ['RFIm-ON-OFF-Noise-mean'];    
+print(hFig, [figDir filesep figFileName '.pdf'], '-dpdf')
+
+%% For the mean traces in fig 3d.
+hFig = createPrintFig(15 * [ 1 2/3]);
+labelsCell = {'Tm4 off x' 'Tm4 off y' 'Tm4 on x' 'Tm4 on y'...
+             'Tm9 off x' 'Tm9 off y' 'Tm9 on x' 'Tm9 on y'};
+rfsToPlot = filteredRFs; % or normRFs.
+plotONOFFTracesFromCell(rfsToPlot, labelsCell, [figDir filesep], brewermap(8, 'Paired'))
+hSubAx = plotLasagnaTracesFromCell(gcf, rfsToPlot, brewermap(8, 'Paired'));
+ylim(hSubAx(:, 1), [-1 1]);
+arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx(1, :))
+cellfun(@(x) caxis(ZeroCenteredBounds(x, [0 100])), rfsToPlot, 'uni', 0)
+arrayfun(@prettifyAxes, hSubAx)
+arrayfun(@offsetAxes, hSubAx)
+set(hSubAx, 'Visible', 'off')
+% figFileName = ['RFcurves-ON-OFF-Noise-XY-lasagna'];    
+% print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+%%
+dataFileName = [figDir filesep 'Fig3.xls'];
+
+DataTable = createRFsTable(rfsToPlot{1, 1}, 'Tm4_OFF_horizontal_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A1')
+DataTable = createRFsTable(rfsToPlot{2, 1}, 'Tm4_OFF_vertical_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A29')
+DataTable = createRFsTable(rfsToPlot{3, 1}, 'Tm4_ON_horizontal_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A56')
+DataTable = createRFsTable(rfsToPlot{4, 1}, 'Tm4_ON_vertical_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A83')
+
+DataTable = createRFsTable(rfsToPlot{1, 2}, 'Tm9_OFF_horizontal_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A110')
+DataTable = createRFsTable(rfsToPlot{2, 2}, 'Tm9_OFF_vertical_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A137')
+DataTable = createRFsTable(rfsToPlot{3, 2}, 'Tm9_ON_horizontal_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A164')
+DataTable = createRFsTable(rfsToPlot{4, 2}, 'Tm9_ON_vertical_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A191')
+
+%%
+splitSign = 1;
+doSmoothing = 1;
+rSquareCutoff = 0.5;
+plotImage = 1;
+hue = [];
+useSubplots = 1;
+% for iChannel = 1: 2
+%     rSquareArray = rsquare{iChannel}';
+%     % Check good fit for all stim.
+%     indsToPlot = all([rSquareArray(:, [1:2]) >= rSquareCutoff ... 
+%                       rSquareArray(:, [3:4]) >= 0.1], 2);
+%     rfInds{:, iChannel} = indsToPlot;
+% end
+
+
+for iChannel = 1: 2
+    rSquareArray = rsquare{iChannel}';
+    % Check good fit for all stim.
+    [rfInds{1: 2, iChannel}] =  deal(all(rSquareArray(:, [1: 2]) >= rSquareCutoff, 2));
+    [rfInds{3: 4, iChannel}] =  deal(all([rSquareArray(:, [1:2]) >= rSquareCutoff ... 
+                                          rSquareArray(:, [3:4]) >= 0.2], 2));
+end
+
+
+%% Figure 3e
+hFig = createPrintFig(15 * [1/2 1/2]);
+
+meanWidth = {mean(width{1}(1: 2, rfInds{1, 1}))
+            mean(width{2}(1: 2, rfInds{1, 2}))
+            mean(width{2}(3: 4, rfInds{3, 2}))};
+
+beanPlot(meanWidth, {'Tm4 off', 'Tm9 off', ' Tm9 on'}, flipud(brewermap(4, 'Set1')), 1, gca,'unbounded', 'vertical', false, true)
+figFileName = ['RF-ON-OFF-Noise-Width-Bean'];    
+print(hFig, [figDir filesep figFileName '.pdf'], '-dpdf')
+
+%% Save data for Fig 3e
+sheetName = 'Fig3e';
+xlswrite(dataFileName, {'Tm4-FHWM-OFF-bars'}, sheetName, 'A1')
+xlswrite(dataFileName, meanWidth{1}, sheetName, 'B1')
+xlswrite(dataFileName, {'Tm9-FHWM-OFF-bars'}, sheetName, 'A2')
+xlswrite(dataFileName, meanWidth{2}, sheetName, 'B2')
+xlswrite(dataFileName, {'Tm9-FHWM-ON-bars'}, sheetName, 'A3')
+xlswrite(dataFileName, meanWidth{3}, sheetName, 'B3')
+
+%% For revision, reviewer 3 plot some raw traces.
+% 
+receptiveFieldTraces = collectRFTraces(dataEdges);
+
+filteredRFTraces = cellfun(@(x, y) x(y, :), receptiveFieldTraces, rfInds, 'uni', 0);
+rng(1001)
+subsampleInds = randsample(min(cellfun(@(x) size(x, 1), receptiveFieldTraces(1, :))) - 1, 10);
+filteredRFTraces = receptiveFieldTraces;
+
+%% Attempt to plot traces correctly, this is working now!
+
+deleteLastROI = 1;
+
+% Define the array to compare to randomized recordings.
+xOffStimFileName = 'StandingStripe_1s_XAxis_5degWide_2degSep_m1.0Con_rand_USEFRUSTUM.txt';
+yOffStimFileName = 'StandingStripe_1s_YAxis_5degWide_2degSep_m1.0Con_rand_USEFRUSTUM.txt';
+xOnStimFileName = 'StandingStripe_1s_XAxis_5degWide_2degSep_p1.0Con_rand_USEFRUSTUM.txt';
+yOnStimFileName = 'StandingStripe_1s_YAxis_5degWide_2degSep_p1.0Con_rand_USEFRUSTUM.txt';
+barStimFileNames = {xOffStimFileName, yOffStimFileName, xOnStimFileName, yOnStimFileName};
+
+%%
+% Take rows as ROIs.
+hFunYOffsetTraces = @(traces, yOffset) ...
+                     bsxfun(@plus, traces, yOffset * (1: size(traces, 1))');
+%%
+nChannels = numel(dataEdges{1}.roiTCsInterp);
+% receptiveFieldsAll = cell(4, nChannels);
+cellCounter = zeros(2, 4);
+hFig = createPrintFig(30 * [1 1 / 2]);
+axInds(1, :) = [1, 2, 5, 6];
+axInds(2, :) = [3, 4, 7, 8];
+for iFly = 1: size(dataEdges, 1)
+%     flyPathInd = iFly;
+    for jStim = 2:5
+        % Find index of matching stimulus.
+        rfInd = find(strcmp(dataEdges{iFly, 1, jStim}.stimulusName, barStimFileNames));
+        if ~isempty(rfInd)
+            for kChannel = 1: nChannels
+                dataTmp = dataEdges{iFly, 1, jStim};
+                nEpochs = numel(dataTmp.stimFrameInds); 
+                nCells = size(dataTmp.roiTCs{kChannel}, 2);
+                roiTCs = dataTmp.roiTCs{kChannel};
+                % Add a fake epoch at the end.
+                dataTmp.stimFrameInds{end + 1} = dataTmp.stimFrameInds{end}(end) + ...
+                                                 (1: numel(dataTmp.stimFrameInds{end}));
+                roiTCs(:, :, dataTmp.stimFrameInds{end}) = nan;                             
+                indsCellArray = [dataTmp.stimFrameInds(1: end - 1); 
+                                 repmat(dataTmp.stimFrameInds(end), 1, nEpochs)];
+                interleavedCellArray = reshape(indsCellArray, 1, []);
+                plotIndsVec = cat(2, interleavedCellArray{:});
+                roiTraces = squeeze(roiTCs(:, 1: end - 1 * deleteLastROI, plotIndsVec));
+                times = dataTmp.frameDuration * (1: numel(plotIndsVec));
+                hAx(kChannel, rfInd) = subplot(2, 4, axInds(kChannel, rfInd));
+                plot(times, hFunYOffsetTraces(roiTraces + cellCounter(kChannel, rfInd), 1));
+                hold on
+                cellCounter(kChannel, rfInd) = cellCounter(kChannel, rfInd) + nCells;
+                rawTracesStruct(iFly, rfInd, kChannel).times = times;
+                rawTracesStruct(iFly, rfInd, kChannel).traces = roiTraces;
+            end
+        end
+    end
+end
+
+%% This is figure S2
+close all
+
+hFig = createPrintFig(30 * [1 1 / 2]);
+axInds(1, :) = [1, 2, 5, 6];
+axInds(2, :) = [3, 4, 7, 8];
+flyInd = 10; % Good examples 4, 6, 7, and 10
+genotypeColors = [166,118,29; 27,158,119] / 255;
+for iStim = 1: 4
+    for jChannel = 1: 2
+        hTracesAx(jChannel, iStim) = subplot(2, 4, axInds(jChannel, iStim));
+        tmpTimes = rawTracesStruct(flyInd, iStim, jChannel).times;
+        tmpTraces = rawTracesStruct(flyInd, iStim, jChannel).traces;
+        plot(tmpTimes, hFunYOffsetTraces(tmpTraces, 0.7), ...
+             'Color', genotypeColors(jChannel, :), 'LineWidth', 1);
+        axis tight
+        hLines = arrayfun(@(x) refline(0, x), 0.7 * (1: size(tmpTraces, 1)));
+        [hLines.Color] = deal([1 1 1] * 0.5);
+        [hLines.LineWidth] = deal(0.5);
+    end
+end
+
+
+xlabel(hTracesAx(5), 'Time (s)')
+ylabel(hTracesAx(5), '\DeltaF/F_0')
+
+setFontForThesis(hTracesAx, gcf);
+arrayfun(@prettifyAxes, hTracesAx)
+arrayfun(@offsetAxes, hTracesAx)
+figFileName = ['Fig3-suppl-meanTraces-sortedPositions'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+%%
+DataTable = createRFsTable(rfsToPlot{1, 1}, 'Tm4_OFF_horizontal_position_deg');
+writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A1')
+%%
+flyInd = 10; % Good examples 4, 6, 7, and 10
+sheetName = 'SupFig3ac';
+startRow = 1;
+stimNameStrCell = {'Tm4_OFF_horizontal_bar_time_seconds', 'Tm4_OFF_vertical_bar_time_seconds', ...
+                   'Tm4_ON_horizontal_bar_time_seconds', 'Tm4_ON_vertical_bar_time_seconds'};
+for iStim = 1: 4
+    DataTable = createRFsTable(rawTracesStruct(flyInd, iStim, 1).traces, stimNameStrCell{iStim});
+    DataTable(:, 1).(stimNameStrCell{iStim}) = rawTracesStruct(flyInd, iStim, 1).times(:);
+    writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
+    startRow = startRow + size(DataTable, 1) + 2;
+end
+
+sheetName = 'SupFig3bd';
+startRow = 1;
+stimNameStrCell = {'Tm9_OFF_horizontal_bar_time_seconds', 'Tm9_OFF_vertical_bar_time_seconds', ...
+                   'Tm9_ON_horizontal_bar_time_seconds', 'Tm9_ON_vertical_bar_time_seconds'};
+for iStim = 1: 4
+    DataTable = createRFsTable(rawTracesStruct(flyInd, iStim, 2).traces, stimNameStrCell{iStim});
+    DataTable(:, 1).(stimNameStrCell{iStim}) = rawTracesStruct(flyInd, iStim, 2).times(:);
+    writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
+    startRow = startRow + size(DataTable, 1) + 2;
+end

+ 200 - 0
figures/MakeTmNeuronsRFs_Fig4.m

@@ -0,0 +1,200 @@
+close all;
+clear all;
+clc;
+
+%% Define search path.
+parentDir = fileparts(fileparts(mfilename('fullpath')));
+saveFolder = [parentDir filesep 'data'];
+load([saveFolder filesep 'Fig4-Tm1Tm2Tm4-processedTableWithBackground'], 'ProcessedTable');
+
+%%
+figDir = [parentDir filesep 'figures' filesep 'Fig4' filesep];
+if ~exist(figDir, 'dir'); mkdir(figDir); end
+
+%%
+% Stimuli of interest.
+stimSetCellStr = getStimSetCellStr(2);
+%% Split table into cell types.
+tokens = arrayfun(@(x)  regexp(x, '.*Tm(\d).*', 'tokens'), ProcessedTable.timeSeriesPath);
+cellTypeNumArray = cellfun(@(x) str2num(x{1}{1}), tokens);
+cellTypeNum = unique(cellTypeNumArray);
+for iNeuron = 1: numel(cellTypeNum)
+    StimTableNeuronCell{iNeuron} = ProcessedTable(cellTypeNumArray == cellTypeNum(iNeuron), :);
+end
+%%
+stimInd = 1;
+nStims = numel(stimSetCellStr);
+for iNeuron = 1: numel(cellTypeNum)
+    for jStim = 1: nStims
+        StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
+    end
+end
+
+for iFly = 1: numel(unique(ProcessedTable.flyInd))
+    flyIndsCell = ProcessedTable.flyInd == iFly;
+    flyStims = ProcessedTable(flyIndsCell, :).stimParamFileName;
+    % Find if stim is within the chosen subset.
+    flyStimInds{iFly} = (cellfun(@(x) find(strcmp(x, stimSetCellStr)), flyStims, 'uni', 0));
+    % If it contains stimuli, check that it includes the 5 chosen.
+    if ~isempty(flyStimInds{iFly})
+        setDiff{iFly} = setdiff(1: 5, [flyStimInds{iFly}{:}]);
+    else
+        setDiff{iFly} = -1;
+    end
+    % For flies with all chosen stimuli order the table accordingly.
+    if isempty(setDiff{iFly})
+        sortFlyIndsTemp = find(flyIndsCell);
+        [~, sortInd] = sort([flyStimInds{iFly}{:}]);
+        sortFlyInds{iFly} = sortFlyIndsTemp(sortInd);
+    end
+end
+
+StimOrderedTable = ProcessedTable(vec([sortFlyInds{:}]), :);
+
+%% Split table into cell types.
+
+tokens = arrayfun(@(x)  regexp(x, '.*Tm(\d).*', 'tokens'), StimOrderedTable.timeSeriesPath);
+cellTypeNumArray = cellfun(@(x) str2num(x{1}{1}), tokens);
+cellTypeNum = unique(cellTypeNumArray);
+for iNeuron = 1: numel(cellTypeNum)
+    StimTableNeuronCell{iNeuron} = StimOrderedTable(cellTypeNumArray == cellTypeNum(iNeuron), :);
+end
+%%
+nStims = numel(stimSetCellStr);
+for iNeuron = 1: numel(cellTypeNum)
+    for jStim = 1: nStims
+        StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
+    end
+end
+
+%% Now can use this table cell array into the other functions.
+
+StimTable = ProcessedTable;
+StimSubTablesCell = arrayfun(@(x) getStimSubsetTable(ProcessedTable, ...
+                                                     stimSetCellStr, x), ...
+                             1: numel(stimSetCellStr), 'uni', 0);  
+
+%% Hand stitched data for Tmneurons and Tm9 saved.
+load([saveFolder filesep 'Fig4-TmNeurons-allNeuronsNeuronStimTableOnOff'], 'tmAll');
+StimTableNeuronStim = tmAll;
+%% Plot tuning curves. Do not discard the on.
+
+alignToOff = true;
+alignToAll = false;
+respQualityThreshold = 0.5;
+rSquareThreshold = 0.2;
+padWithNaNs = 0;
+qualityMeasure = 'responseQuality';
+
+%%
+cellTypeStr = {'Tm1', 'Tm2', 'Tm4', 'Tm9'};
+% plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), 0.0, 0.5, figDir, [cellTypeStr{:}]);
+plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), 0.2, 0.5, figDir, [cellTypeStr{:}]);
+% plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), -20.0, 0, figDir, [cellTypeStr{:}]);
+%% Export excel data figure 4b-i
+receptiveFields = plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), 0.2, 0.5, figDir, [cellTypeStr{:}]);
+%% Write data per cell type
+dataFileName = [figDir filesep 'Fig4.xls'];
+sheetNamesCell = {'bf', 'cg', 'dh', 'ei'};
+
+stimNameStrCell = {'_OFF_horizontal_position_deg', '_OFF_vertical_position_deg', ...
+                   '_ON_horizontal_position_deg', '_ON_vertical_position_deg'};
+               
+for iCellType = 1: numel(cellTypeStr)
+    startRow = 1;
+    sheetName = ['Fig4' sheetNamesCell{iCellType}];
+    for iStim = 1: 4
+        parameterLabel = [cellTypeStr{iCellType} stimNameStrCell{iStim}];
+        DataTable = createRFsTable(receptiveFields{iCellType, iStim}, parameterLabel);
+        writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
+        startRow = startRow + size(DataTable, 1) + 2;
+    end
+end
+
+%% Export data from Supplementary Figure 4
+receptiveFields = plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), -20, 0, figDir, [cellTypeStr{:}]);
+%%
+dataFileName = [figDir filesep 'Fig4.xls'];
+sheetNamesCell = {'ae', 'bf', 'cg', 'dh'};
+
+stimNameStrCell = {'_OFF_horizontal_position_deg', '_OFF_vertical_position_deg', ...
+                   '_ON_horizontal_position_deg', '_ON_vertical_position_deg'};
+               
+for iCellType = 1: numel(cellTypeStr)
+    startRow = 1;
+    sheetName = ['SuppFig4' sheetNamesCell{iCellType}];
+    for iStim = 1: 4
+        parameterLabel = [cellTypeStr{iCellType} stimNameStrCell{iStim}];
+        DataTable = createRFsTable(receptiveFields{iCellType, iStim}, parameterLabel);
+        writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
+        startRow = startRow + size(DataTable, 1) + 2;
+    end
+end
+
+%% Inline from plotRFsAlignedWithFitGenotypes
+respIndThresh = respQualityThreshold;
+rSqThresh = rSquareThreshold;
+clear flyIndsCell
+for iNeuron = 1: size(StimTableNeuronStim(:, 2: end), 2)
+    StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
+    [~, responseIndArrayCell, rSquareArrayCell, fitParamsCell, ...
+      flyIndsCell{iNeuron}, ~, ~] = getStimArraysForSameRois(StimTableCell);                                          
+    % Align to fit
+    widthCell{iNeuron} = cellfun(@(x) arrayfun(@(y) 2 * (2 * sqrt(log(2))) * y.c1, x, 'uni', 0), ...
+                        fitParamsCell, 'uni', 1); % The 2 factor is because of bar separation.
+    validIndsCell{iNeuron} = cellfun(@(x, y) x(:) > rSqThresh & y(:) > respIndThresh, ...
+                                     rSquareArrayCell, responseIndArrayCell, 'uni', 0)';
+end
+%%
+validOffInds = cellfun(@(x) x{1} & x{2}, validIndsCell, 'uni', 0);
+validOnInds = cellfun(@(x)  x{3} & x{4}, validIndsCell, 'uni', 0);
+%%
+widthCell = cellfun(@(x) cat(1, x{:}), widthCell, 'uni', 0);
+%%
+xOFFWidth = cellfun(@(x, y) x(1, y), widthCell, validOffInds, 'uni', 0);
+yOFFWidth = cellfun(@(x, y) x(2, y), widthCell, validOffInds, 'uni', 0);
+xONWidth = cellfun(@(x, y) x(3, y), widthCell, validOnInds, 'uni', 0);
+yONWidth = cellfun(@(x, y) x(4, y), widthCell, validOnInds, 'uni', 0);
+allWidthsCell = {xOFFWidth, yOFFWidth, xONWidth, yONWidth};
+%%
+for iW = 1: numel(allWidthsCell)
+    allWidthsCell{iW}(cellfun(@isempty, allWidthsCell{iW}, 'uni', 1)) = {pi};
+end
+%%
+hFig = createPrintFig(15 * [ 1 1/4]);
+colors = brewermap(4, 'Set1');
+colors = colors([1, 2, 4, 3], :);
+for iAx = 1: numel(allWidthsCell)
+    hSubAx(iAx) = subplot(1, 4, iAx);
+    hWidthAx = beanPlot(allWidthsCell{iAx}, cellTypeStr, colors, ...
+                        1, gca, [0-eps 100], 'vertical', 0, 1, 0);
+end
+linkaxes(hSubAx, 'y');
+setFontForThesis(hWidthAx, gcf);
+arrayfun(@prettifyAxes, hWidthAx)
+arrayfun(@offsetAxes, hWidthAx)
+figFileName = ['RF-FWHM-bean'];    
+print(hFig, [figDir figFileName '.pdf'], '-dpdf')
+%% Export data from figure 4 j-m
+%%
+dataFileName = [figDir filesep 'Fig4.xls'];
+
+stimNameStrCell = {'_OFF_bars_horizontal_deg', '_OFF_bars_vertical_deg', ...
+                   '_ON_bars_horizontal_deg', '_ON_bars_vertical_deg'};
+startRow = 1;
+sheetName = ['Fig4j-m'];
+for iCellType = 1: numel(cellTypeStr)
+    for iStim = 1: 4
+        parameterLabel = [cellTypeStr{iCellType} '_FWHM' stimNameStrCell{iStim}];
+        xlswrite(dataFileName, {parameterLabel}, sheetName, ['A' num2str(startRow)])
+        if allWidthsCell{iStim}{iCellType} == pi % pi was placeholder for no data
+            xlswrite(dataFileName, {''}, sheetName, ['B' num2str(startRow)])
+        else
+            xlswrite(dataFileName, allWidthsCell{iStim}{iCellType}, sheetName, ['B' num2str(startRow)])
+        end
+        startRow = startRow + 1;
+    end
+    startRow = startRow + 1;
+end
+%%
+meanONWidth = cellfun(@mean, cellfun(@(x, y) (x + y) / 2, xONWidth, yONWidth, 'uni', 0));