123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719 |
- 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
|