MakeFigTmNeuronsCDMTableRevisionDoGFitConstraint_Fig7.m 46 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103
  1. close all;
  2. clear all;
  3. clc;
  4. %% Define search path.
  5. parentDir = fileparts(fileparts(mfilename('fullpath')));
  6. saveFolder = [parentDir filesep 'data'];
  7. load([saveFolder filesep 'Fig7-Tm9CDM-processedTableWithoutBackgroundOnOff'], 'ProcessedTable');
  8. %%
  9. figDir = [parentDir filesep 'figures' filesep 'Fig7' filesep];
  10. if ~exist(figDir, 'dir'); mkdir(figDir); end
  11. %%
  12. % Stimuli of interest.
  13. stimSetCellStr = getStimSetCellStr(2);
  14. % Check only bars.
  15. nStims = 5;
  16. stimSetCellStr = stimSetCellStr(1: nStims);
  17. %% Split table into cell types.
  18. tokens = arrayfun(@(x) regexp(x, '.*Tm(\d).*', 'tokens'), ProcessedTable.timeSeriesPath);
  19. cellTypeNumArray = cellfun(@(x) str2num(x{1}{1}), tokens);
  20. cellTypeNum = unique(cellTypeNumArray);
  21. for iNeuron = 1: numel(cellTypeNum)
  22. StimTableNeuronCell{iNeuron} = ProcessedTable(cellTypeNumArray == cellTypeNum(iNeuron), :);
  23. end
  24. %%
  25. stimInd = 1;
  26. nStims = numel(stimSetCellStr);
  27. for iNeuron = 1: numel(cellTypeNum)
  28. for jStim = 1: nStims
  29. StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
  30. end
  31. end
  32. for iFly = 1: numel(unique(ProcessedTable.flyInd))
  33. flyInds = ProcessedTable.flyInd == iFly;
  34. flyStims = ProcessedTable(flyInds, :).stimParamFileName;
  35. % Find if stim is within the chosen subset.
  36. flyStimInds{iFly} = (cellfun(@(x) find(strcmp(x, stimSetCellStr)), flyStims, 'uni', 0));
  37. % If it contains stimuli, check that it includes the 5 chosen.
  38. if ~isempty(flyStimInds{iFly})
  39. setDiff{iFly} = setdiff(1: nStims, [flyStimInds{iFly}{:}]);
  40. else
  41. setDiff{iFly} = -1;
  42. end
  43. % For flies with all chosen stimuli order the table accordingly.
  44. if isempty(setDiff{iFly})
  45. sortFlyIndsTemp = find(flyInds);
  46. [~, sortInd] = sort([flyStimInds{iFly}{:}]);
  47. sortFlyInds{iFly} = sortFlyIndsTemp(sortInd);
  48. end
  49. end
  50. StimOrderedTable = ProcessedTable(vec(cat(1, sortFlyInds{:})), :);
  51. %%
  52. %% Split table into cell types.
  53. tokens = arrayfun(@(x) regexp(x, '.*Tm(\d).*plus(.*CDM)', 'tokens'), StimOrderedTable.timeSeriesPath);
  54. tokens = cellfun(@(x) [x{1}{:}], tokens, 'uni', 0);
  55. cellTypeNum = unique(tokens);
  56. for iNeuron = 1: numel(cellTypeNum)
  57. StimTableNeuronCell{iNeuron} = StimOrderedTable(strcmp(tokens, cellTypeNum(iNeuron)), :);
  58. end
  59. %%
  60. nStims = numel(stimSetCellStr);
  61. for iNeuron = 1: numel(cellTypeNum)
  62. for jStim = 1: nStims
  63. StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
  64. end
  65. end
  66. %%
  67. % Maybe this should get its own file
  68. doFit = 0;
  69. if doFit
  70. flyFits = fitDifferenceOfGaussiansFromNeuronStimTables(StimTableNeuronStim);
  71. save([dataSaveFolder 'DoGFitsTm9TestFminCon'], 'flyFits', '-v7.3');
  72. end
  73. %%
  74. load([saveFolder filesep 'Fig7-Tm9CDM-DoGFitsTm9TestFminCon'], 'flyFits');
  75. %% From DoG
  76. if size(flyFits, 2) == 5
  77. flyFits(:, 1, :) = [];
  78. end
  79. cdm = squeeze(flyFits(1, :, :));
  80. noCdm = squeeze(flyFits(2, :, :));
  81. cdm(:, all(cellfun(@isempty, cdm))) = [];
  82. noCdm(:, all(cellfun(@isempty, noCdm))) = [];
  83. %% Bad chunk of code to try to align data structure of main analysis with the
  84. % new Diff of Gaussians fit structure.
  85. rSqThresh = 0;
  86. respIndThresh = 0.5;
  87. iInd = 1;
  88. clear roisToRemoveInd
  89. for iNeuron = [2, 1] % first control then cdm
  90. StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
  91. [tcArrayCell, responseIndArrayCell, rSquareArrayCell, fitParamsCell, ...
  92. flyIndsCell, tcExtCellArray, tcArgExtCellArray, roisToRemoveIndTmp] = ...
  93. getStimArraysForSameRois(StimTableCell, 1);
  94. roisToRemoveInd{iInd} = roisToRemoveIndTmp;
  95. validIndsCell{iInd} = cellfun(@(x, y) x(:) > rSqThresh & y(:) > respIndThresh, ...
  96. rSquareArrayCell, responseIndArrayCell, 'uni', 0)';
  97. iInd = iInd + 1;
  98. end
  99. [flyDelInd, stimDelInd] = cellfun(@(x) ind2sub(size(x), ...
  100. find(~cellfun(@isempty, x))), ...
  101. roisToRemoveInd, 'uni', 0);
  102. %% for control
  103. genotypeInd = 1;
  104. roiDelInds = roisToRemoveInd{genotypeInd};
  105. delIndsPerFly = {};
  106. for iRoi = 1: numel(flyDelInd{genotypeInd})
  107. delInds = roiDelInds{flyDelInd{genotypeInd}(iRoi), ...
  108. stimDelInd{genotypeInd}(iRoi)}{1};
  109. delIndsPerFly{flyDelInd{genotypeInd}(iRoi)} = [];
  110. delIndsPerFly{flyDelInd{genotypeInd}(iRoi)} = [delIndsPerFly{flyDelInd{genotypeInd}(iRoi)}, delInds];
  111. end
  112. for iFly = 1: numel(delIndsPerFly)
  113. if ~isempty(delIndsPerFly{iFly})
  114. for jStim = 1: size(noCdm, 1)
  115. tmp = noCdm{jStim, iFly};
  116. tmp(unique(delIndsPerFly{iFly})) = [];
  117. noCdm{jStim, iFly} = tmp;
  118. end
  119. end
  120. end
  121. %% for cdm
  122. clear delIndsPerFly
  123. genotypeInd = 2;
  124. roiDelInds = roisToRemoveInd{genotypeInd};
  125. delIndsPerFly = {};
  126. for iRoi = 1: numel(flyDelInd{genotypeInd})
  127. delInds = roiDelInds{flyDelInd{genotypeInd}(iRoi), ...
  128. stimDelInd{genotypeInd}(iRoi)}{1};
  129. delIndsPerFly{flyDelInd{genotypeInd}(iRoi)} = [];
  130. delIndsPerFly{flyDelInd{genotypeInd}(iRoi)} = [delIndsPerFly{flyDelInd{genotypeInd}(iRoi)}, delInds];
  131. end
  132. for iFly = 1: numel(delIndsPerFly)
  133. if ~isempty(delIndsPerFly{iFly})
  134. for jStim = 1: size(cdm, 1)
  135. tmp = cdm{jStim, iFly};
  136. tmp(unique(delIndsPerFly{iFly})) = [];
  137. cdm{jStim, iFly} = tmp;
  138. end
  139. end
  140. end
  141. %%
  142. % reviewer #1 analysis by fly
  143. % this method is not aligned with the new fit structure, but the structure
  144. % is separated by flies already so we can recreate the indices.
  145. % the minux is to exclude the background roi.
  146. flyIndsCdm = repelem(1: size(cdm, 2), cellfun(@numel, cdm(1, :)) - 1);
  147. flyIndsNoCdm = repelem(1: size(noCdm, 2), cellfun(@numel, noCdm(1, :)) - 1);
  148. %%
  149. [Amp, Pos, Sigma, RSq, baseline] = deal(cell(1, 2));
  150. [Amp{1}, Pos{1}, Sigma{1}, RSq{1}, baseline{1}] = parseFitParams(noCdm);
  151. [Amp{2}, Pos{2}, Sigma{2}, RSq{2}, baseline{2}] = parseFitParams(cdm);
  152. centerSurround = cellfun(@(x, y) x ./ y, Amp{1}.center, Amp{1}.surround, 'uni', 0);
  153. %%
  154. allAmpCenter = cellfun(@(x) cell2mat(x.center), Amp, 'uni', 0);
  155. allAmpSurround = cellfun(@(x) cell2mat(x.surround), Amp, 'uni', 0);
  156. allAmpCenterSurround = cellfun(@(x, y) x ./ y, allAmpCenter, allAmpSurround, 'uni', 0);
  157. if any(~cellfun(@isnan, baseline))
  158. allAmpBaseline = cellfun(@(x) cell2mat(x), baseline, 'uni', 0);
  159. end
  160. allRSq = cellfun(@cell2mat, RSq, 'uni', 0);
  161. %%
  162. minOffRSq = cellfun(@(x) min(x(1:2, :)), allRSq, 'uni', 0);
  163. minOnRSq = cellfun(@(x) min(x(3:4, :)), allRSq, 'uni', 0);
  164. minAllRSq = cellfun(@min, allRSq, 'uni', 0);
  165. %%
  166. respIndThresh = 0.5;
  167. % thresholdRSq = 0.2;
  168. thresholdRSq = -20;
  169. goodCells = cellfun(@(x) x > thresholdRSq, minAllRSq, 'uni', 0);
  170. %% This is to plot matching RFs to dog fit
  171. % plotRFsAlignedWithFitGenotypes(StimTableNeuronStim([2, 1], 2: end), 0, 0.5, figDir, [cellTypeStr{[2, 1]}]);
  172. [offX, offY, onX, onY] = deal(cell(size(StimTableNeuronStim, 1), 1));
  173. validIndsCellTmp = goodCells([2, 1]); % flip from (noCDM, CDM) to conserve original order in table (CDM, noCDM).
  174. for iNeuron = [2, 1] % 2: noCDM, 1: CDM, thus offX = [cdm, noCDM]
  175. StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
  176. % Overwrite the gauss1 rsquare with gauss2 rsquare.
  177. [tcArrayCell, responseIndArrayCell, ~, fitParamsCell, ...
  178. flyIndsCell, tcExtCellArray, tcArgExtCellArray] = getStimArraysForSameRois(StimTableCell);
  179. % Align to fit
  180. posCell = cellfun(@(x) arrayfun(@(y) y.b1, x, 'uni', 0), fitParamsCell, 'uni', 1);
  181. lagsCell = cellfun(@(x, y) round(x - size(y, 2) / 2), posCell, tcArrayCell, 'uni', 0);
  182. barRFs = cellfun(@(x, y) alignRFsToAbsMax(x, 0, y'), tcArrayCell, lagsCell, 'uni', 0);
  183. %% Blank epoch was already deleted, now delete vertical bars on edges of screen.
  184. barRFs{2}(:, end-1: end) = [];
  185. barRFs{4}(:, end-1: end) = [];
  186. barRFs{2}(:, 1: 2) = [];
  187. barRFs{4}(:, 1: 2) = [];
  188. %%
  189. % Include the response quality index.
  190. validRespInds = all(cat(2, responseIndArrayCell{:}) > respIndThresh, 2)';
  191. offX{iNeuron} = barRFs{1}(validIndsCellTmp{iNeuron} & validRespInds, :);
  192. offY{iNeuron} = barRFs{2}(validIndsCellTmp{iNeuron} & validRespInds, :);
  193. onX{iNeuron} = barRFs{3}(validIndsCellTmp{iNeuron} & validRespInds, :);
  194. onY{iNeuron} = barRFs{4}(validIndsCellTmp{iNeuron} & validRespInds, :);
  195. validRespIndsCell{iNeuron} = validRespInds;
  196. end
  197. % Include the response quality index
  198. goodCells = cellfun(@(x, y) x & y, goodCells, validRespIndsCell([2, 1]), 'uni', 0);
  199. %%
  200. [offX{cellfun(@isempty, offX)}] = deal(zeros(1, size(offX{find(cellfun(@isempty, offX)==0, 1)}, 2)));
  201. [offY{cellfun(@isempty, offY)}] = deal(zeros(1, size(offY{find(cellfun(@isempty, offY)==0, 1)}, 2)));
  202. [onX{cellfun(@isempty, onX)}] = deal(zeros(1, size(onX{find(cellfun(@isempty, onX)==0, 1)}, 2)));
  203. [onY{cellfun(@isempty, onY)}] = deal(zeros(1, size(onY{find(cellfun(@isempty, onY)==0, 1)}, 2)));
  204. %% Plot selected rf traces for figure 7a-b
  205. cellTypeStr = {'CDM', 'noCDM'};
  206. plotLasagnaRFs(offX([2, 1]), offY([2, 1]), 'OFF', figDir, [cellTypeStr{[2, 1]}]);
  207. plotLasagnaRFs(onX([2, 1]), onY([2, 1]), 'ON', figDir, [cellTypeStr{[2, 1]}]);
  208. %% Export data from figure 7a-b
  209. dataFileName = [figDir filesep 'Fig7.xls'];
  210. stimNameStrCell = {'_OFF_horizontal_position_deg', '_OFF_vertical_position_deg', ...
  211. '_ON_horizontal_position_deg', '_ON_vertical_position_deg'};
  212. cellTypeStr = {'CDM', 'Control'};
  213. receptiveFields = [offX(1) offY(1) onX(1) onY(1); offX(2) offY(2) onX(2) onY(2)];
  214. startRow = 1;
  215. sheetName = ['Fig7ab'];
  216. for iCellType = 1: numel(cellTypeStr)
  217. for iStim = 1: 4
  218. parameterLabel = [cellTypeStr{iCellType} stimNameStrCell{iStim}];
  219. DataTable = createRFsTable(receptiveFields{iCellType, iStim}, parameterLabel);
  220. writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
  221. startRow = startRow + size(DataTable, 1) + 2;
  222. end
  223. end
  224. %% For reviewer #1 comment 19. Plotting traces with different normalizations.
  225. % this normalizes each tuning curve, so width differences are highlighted
  226. % regardless of amplitude differences between conditions.
  227. % normalize to peak and smooth and normalize again.
  228. maxAbsNorm = @(x, dim) x ./ max(abs(x), [], dim);
  229. maxNorm = @(x, dim) x ./ max(x, [], dim);
  230. minNorm = @(x, dim) x ./ -min(x, [], dim);
  231. rangeNorm = @(x, dim) x ./ range(x, dim);
  232. maxNormRFs = cellfun(@(x) maxNorm(x, 2), [offX, offY], 'uni', 0);
  233. minNormRFs = cellfun(@(x) minNorm(x, 2), [onX, onY], 'uni', 0);
  234. [normOffX, normOffY, normOnX, normOnY] = deal(maxNormRFs([2, 1], 1), maxNormRFs([2, 1], 2), ...
  235. minNormRFs([2, 1], 1), minNormRFs([2, 1], 2));
  236. plotLasagnaRFs(normOffX, normOffY, 'OFF', figDir, [cellTypeStr{[2, 1]} '-maxMinNorm']);
  237. plotLasagnaRFs(normOnX, normOnY, 'ON', figDir, [cellTypeStr{[2, 1]} '-maxMinNorm']);
  238. rangeNormRFs = cellfun(@(x) rangeNorm(x, 2), [offX, offY, onX, onY], 'uni', 0);
  239. [normOffX, normOffY, normOnX, normOnY] = deal(rangeNormRFs([2, 1], 1), rangeNormRFs([2, 1], 2), ...
  240. rangeNormRFs([2, 1], 3), rangeNormRFs([2, 1], 4));
  241. plotLasagnaRFs(normOffX, normOffY, 'OFF', figDir, [cellTypeStr{[2, 1]} '-rangeNorm']);
  242. plotLasagnaRFs(normOnX, normOnY, 'ON', figDir, [cellTypeStr{[2, 1]} '-rangeNorm']);
  243. %% Actually quantification is averaging over orientations, so let's plot those
  244. normOff = cellfun(@(x, y) (x + y) / 2, normOffX, normOffY, 'uni', 0);
  245. normOn = cellfun(@(x, y) (x + y) / 2, normOnX, normOnY, 'uni', 0);
  246. plotLasagnaRFs(normOff, normOn, 'ONOFF', figDir, [cellTypeStr{[2, 1]} '-rangeNorm-meanXY']);
  247. %%
  248. offAmpCenter = cellfun(@(x, y) mean(x(1:2, y)), allAmpCenter, goodCells, 'uni', 0);
  249. offAmpSurround = cellfun(@(x, y) mean(x(1:2, y)), allAmpSurround, goodCells, 'uni', 0);
  250. onAmpCenter = cellfun(@(x, y) mean(x(3:4, y)), allAmpCenter, goodCells, 'uni', 0);
  251. onAmpSurround = cellfun(@(x, y) mean(x(3:4, y)), allAmpSurround, goodCells, 'uni', 0);
  252. %%
  253. sigma2fwhmFactor = 2 * 2 * sqrt(log(2));
  254. allWidthCenter = cellfun(@(x) cell2mat(x.center) * sigma2fwhmFactor, Sigma, 'uni', 0);
  255. allWidthSurround = cellfun(@(x) cell2mat(x.surround) * sigma2fwhmFactor, Sigma, 'uni', 0);
  256. % allWidthCenterSurround = cellfun(@(x, y) x ./ y, allWidthCenter, allWidthSurround, 'uni', 0);
  257. %%
  258. offWidthCenter = cellfun(@(x, y) mean(x(1:2, y)), allWidthCenter, goodCells, 'uni', 0);
  259. offWidthSurround = cellfun(@(x, y) mean(x(1:2, y)), allWidthSurround, goodCells, 'uni', 0);
  260. onWidthCenter = cellfun(@(x, y) mean(x(3:4, y)), allWidthCenter, goodCells, 'uni', 0);
  261. onWidthSurround = cellfun(@(x, y) mean(x(3:4, y)), allWidthSurround, goodCells, 'uni', 0);
  262. %% Start plotting
  263. flyIndsCell = {flyIndsNoCdm, flyIndsCdm};
  264. goodFlyIndsCell = cellfun(@(x, y) x(y), flyIndsCell, goodCells, 'uni', 0);
  265. offWidthCenterFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
  266. flyIndsCell, offWidthCenter, goodCells, ...
  267. 'uni', 0);
  268. offWidthSurroundFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
  269. flyIndsCell, offWidthSurround, goodCells, ...
  270. 'uni', 0);
  271. onWidthCenterFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
  272. flyIndsCell, onWidthCenter, goodCells, ...
  273. 'uni', 0);
  274. onWidthSurroundFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
  275. flyIndsCell, onWidthSurround, goodCells, ...
  276. 'uni', 0);
  277. offAmpCenterFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
  278. flyIndsCell, offAmpCenter, goodCells, ...
  279. 'uni', 0);
  280. offAmpSurroundFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
  281. flyIndsCell, offAmpSurround, goodCells, ...
  282. 'uni', 0);
  283. onAmpCenterFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
  284. flyIndsCell, onAmpCenter, goodCells, ...
  285. 'uni', 0);
  286. onAmpSurroundFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
  287. flyIndsCell, onAmpSurround, goodCells, ...
  288. 'uni', 0);
  289. emptyFlyInds = cellfun(@(x, y) find(accumarray(x(y)', 1) == 0), flyIndsCell, goodCells, 'uni', 0);
  290. % Remove empty flies form array for plotting and quantification.
  291. for iNeuron = 1: 2
  292. onAmpCenterFlyMeans{iNeuron}(emptyFlyInds{iNeuron}) = [];
  293. offAmpCenterFlyMeans{iNeuron}(emptyFlyInds{iNeuron}) = [];
  294. onAmpSurroundFlyMeans{iNeuron}(emptyFlyInds{iNeuron}) = [];
  295. offAmpSurroundFlyMeans{iNeuron}(emptyFlyInds{iNeuron}) = [];
  296. end
  297. %% Plot amplitudes of center vs surround. Figure 7e
  298. clear hSubAx
  299. hFig = createPrintFig(10 * [3/2 1/2]);
  300. % suptitle('Receptive field amplitude of center vs surround')
  301. hSubAx(1) = subplot(1, 3, 1);
  302. scatter(offAmpCenter{1}, offAmpSurround{1}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
  303. hold on
  304. scatter(offAmpCenter{2}, offAmpSurround{2}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
  305. lsline
  306. % refline(-1, 0)
  307. axis equal square
  308. title('OFF RF Amplitude')
  309. xlabel('Center (\Delta F / F_0)')
  310. ylabel('Surround (\Delta F / F_0)')
  311. legend({'Control', 'CDM'}, 'Location', 'southwest')
  312. hSubAx(2) = subplot(1, 3, 2);
  313. scatter(onAmpCenter{1}, onAmpSurround{1}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
  314. hold on
  315. scatter(onAmpCenter{2}, onAmpSurround{2}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
  316. lsline
  317. % refline(-1, 0)
  318. axis equal square
  319. title('ON RF Amplitude')
  320. hSubAx(3) = subplot(1, 3, 3);
  321. scatter(onAmpCenter{1}, onAmpSurround{1}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
  322. hold on
  323. scatter(onAmpCenter{2}, onAmpSurround{2}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
  324. lsline
  325. % refline(-1, 0)
  326. axis equal square
  327. xlim([-0.7 -0.2]);
  328. ylim([0 0.5]);
  329. title('ON RF Amplitude')
  330. arrayfun(@prettifyAxes, hSubAx);
  331. setFontForThesis(hSubAx, hFig);
  332. % figDir = 'P:\Documents\Figures\Paper\Tm9\CDM\';
  333. figFileName = ['CenterSurround-Amp-Scatter'];
  334. set(gcf,'renderer','painters')
  335. % set(gcf, 'PaperPositionMode', 'auto');
  336. print(hFig, [figDir figFileName '.pdf'], '-dpdf')
  337. %% Correlations for figure 7e
  338. [offAmpCorr, offAmpCorrP] = cellfun(@(x, y) corr(x(:), y(:)), ...
  339. offAmpCenter, offAmpSurround, 'uni', 1);
  340. [onAmpCorr, onAmpCorrP] = cellfun(@(x, y) corr(x(:), y(:)), ...
  341. onAmpCenter, onAmpSurround, 'uni', 1);
  342. % Bootstrap CI
  343. [offAmpCorrCI, offAmpCorrs] = cellfun(@(x, y) bootci(1000, @corr, x(:), y(:)), ...
  344. offAmpCenter, offAmpSurround, 'uni', 0);
  345. [onAmpCorrCI, onAmpCorrs] = cellfun(@(x, y) bootci(1000, @corr, x(:), y(:)), ...
  346. onAmpCenter, onAmpSurround, 'uni', 0);
  347. offAmpCorrCIArray = cat(2, offAmpCorrCI{:}) % each column is one condition
  348. offAmpCorrsArray = cellfun(@mean, offAmpCorrs)
  349. onAmpCorrCIArray = cat(2, onAmpCorrCI{:}) % each column is one condition
  350. onAmpCorrsArray = cellfun(@mean, onAmpCorrs)
  351. %% Export data for figure 7e.
  352. sheetName = ['Fig7e'];
  353. startRow = 1;
  354. conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
  355. varNames = {[conditionStrCell{1} '_OFF_RF_Center_Amplitude_dF_F0'], ...
  356. [conditionStrCell{1} '_OFF_RF_Surround_Amplitude_dF_F0'], ...
  357. [conditionStrCell{2} '_OFF_RF_Center_Amplitude_dF_F0'], ...
  358. [conditionStrCell{2} '_OFF_RF_Surround_Amplitude_dF_F0']};
  359. DataTable = table(offAmpCenter{1}(:), offAmpSurround{1}(:), ...
  360. 'VariableNames', varNames(1: 2));
  361. writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
  362. DataTable = table(offAmpCenter{2}(:), offAmpSurround{2}(:), ...
  363. 'VariableNames', varNames(3: 4));
  364. writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['D' num2str(startRow)]);
  365. startRow = max(cellfun(@numel, offAmpCenter)) + 3;
  366. conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
  367. varNames = {[conditionStrCell{1} '_ON_RF_Center_Amplitude_dF_F0'], ...
  368. [conditionStrCell{1} '_ON_RF_Surround_Amplitude_dF_F0'], ...
  369. [conditionStrCell{2} '_ON_RF_Center_Amplitude_dF_F0'], ...
  370. [conditionStrCell{2} '_ON_RF_Surround_Amplitude_dF_F0']};
  371. DataTable = table(onAmpCenter{1}(:), onAmpSurround{1}(:), ...
  372. 'VariableNames', varNames(1: 2));
  373. writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
  374. DataTable = table(onAmpCenter{2}(:), onAmpSurround{2}(:), ...
  375. 'VariableNames', varNames(3: 4));
  376. writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['D' num2str(startRow)]);
  377. %% Same but per fly for revision
  378. hFig = createPrintFig(10 * [2/2 1/2]);
  379. % suptitle('Receptive field amplitude of center vs surround')
  380. hSubAx(1) = subplot(1, 2, 1);
  381. scatter(offAmpCenterFlyMeans{1}, offAmpSurroundFlyMeans{1}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
  382. hold on
  383. scatter(offAmpCenterFlyMeans{2}, offAmpSurroundFlyMeans{2}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
  384. lsline
  385. % refline(-1, 0)
  386. axis equal square
  387. title('OFF RF Amplitude')
  388. xlabel('Center (\Delta F / F_0)')
  389. ylabel('Surround (\Delta F / F_0)')
  390. legend({'Control', 'CDM'}, 'Location', 'southwest')
  391. hSubAx(2) = subplot(1, 2, 2);
  392. scatter(onAmpCenterFlyMeans{1}, onAmpSurroundFlyMeans{1}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
  393. hold on
  394. scatter(onAmpCenterFlyMeans{2}, onAmpSurroundFlyMeans{2}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
  395. lsline
  396. % refline(-1, 0)
  397. axis equal square
  398. title('ON RF Amplitude')
  399. arrayfun(@prettifyAxes, hSubAx);
  400. setFontForThesis(hSubAx, hFig);
  401. % figDir = 'P:\Documents\Figures\Paper\Tm9\CDM\';
  402. figFileName = ['CenterSurround-Amp-Scatter-FlyMeans'];
  403. set(gcf,'renderer','painters')
  404. % set(gcf, 'PaperPositionMode', 'auto');
  405. print(hFig, [figDir figFileName '.pdf'], '-dpdf')
  406. %% RF params stats from DoG fits
  407. paramsToTest = {offAmpCenter, offAmpSurround, onAmpCenter, onAmpSurround, ...
  408. offWidthCenter, offWidthSurround, onWidthCenter, onWidthSurround};
  409. %%
  410. nROIsPerFlyCell = cellfun(@(x) accumarray(x(:), 1), goodFlyIndsCell, 'uni', 0);
  411. % Delete flies with no rois to analyze.
  412. nROIsPerFlyCell{1}(nROIsPerFlyCell{1} == 0) = [];
  413. nROIsPerFlyCell{2}(nROIsPerFlyCell{2} == 0) = [];
  414. % generate continuos numbering.
  415. reorderedFlyIndsCell = cellfun(@(x) repelem(1: numel(x), x), nROIsPerFlyCell, 'uni', 0);
  416. % samples are cells, but we will permute all cells of a fly together.
  417. sample1 = 1: numel(nROIsPerFlyCell{1});
  418. sample2 = numel(sample1) + (1: numel(nROIsPerFlyCell{2}));
  419. %%
  420. permutations = 10000;
  421. sidedness = 'both';
  422. plotresult = 1;
  423. pVals = nan(1, numel(paramsToTest));
  424. obsDiffs = nan(1, numel(paramsToTest));
  425. for iParam = 1: numel(paramsToTest)
  426. % Permutation test needs row vectors
  427. paramsQuality = cellfun(@(x) x(:)', paramsToTest{iParam}, 'uni', 0);
  428. [obsDiffs(iParam), pVals(iParam)] = nestetPermutationTest(sample1, sample2, ...
  429. paramsQuality, reorderedFlyIndsCell, ...
  430. permutations, sidedness, plotresult);
  431. end
  432. %% Plotting RF params for figure 7c-d
  433. positiveSupport = [0-eps 10000];
  434. negativeSupport = [-1000 0+eps];
  435. hFig = createPrintFig(15 * [1/2 1/3]);
  436. colors = getONOFFXYColors;
  437. dark2 = brewermap(8, 'Dark2');
  438. colors = dark2([3, 2], :);
  439. % stimCellStr = {'XOFF', 'YOFF', 'XON', 'YON'};
  440. stimCellStr = {'Control', 'CDM'};
  441. titleStrCell = {'OFF Center', 'Surround', 'ON Center', 'Surround'};
  442. dataCell = {offAmpCenter, offAmpSurround, onAmpCenter, onAmpSurround};
  443. supportCell = {positiveSupport, negativeSupport, negativeSupport, positiveSupport};
  444. for iAx = 1: 4
  445. hSubAx(iAx) = subplot(1, 4, iAx);
  446. title(titleStrCell{iAx})
  447. beanPlot(dataCell{iAx}, stimCellStr, colors, 1, gca, supportCell{iAx}, 'vertical', 0);
  448. text(mean(xlim), 0.9 * getArrayElem(ylim, 2), ...
  449. sprintf('p = %1.1e', pVals(iAx)), 'HorizontalAlignment', 'center');
  450. axis tight
  451. end
  452. % linkaxes(hSubAx(1:2), 'xy');
  453. % linkaxes(hSubAx(3:4), 'xy');
  454. ylabel(hSubAx(1), 'RF Amplitude (\Delta F/F_0)')
  455. hFig.Color = 'none';
  456. set(hSubAx, 'Color', 'none');
  457. arrayfun(@prettifyAxes, hSubAx)
  458. setFontForThesis(hSubAx, hFig);
  459. figFileName = ['CenterSurround-Amp-Bean'];
  460. set(gcf,'renderer','painters')
  461. print(hFig, [figDir figFileName '.pdf'], '-dpdf')
  462. hFig = createPrintFig(15 * [1/2 1/3]);
  463. dataCell = {offWidthCenter, offWidthSurround, onWidthCenter, onWidthSurround};
  464. for iAx = 1: 4
  465. hSubAx(iAx) = subplot(1, 4, iAx);
  466. title(titleStrCell{iAx})
  467. beanPlot(dataCell{iAx}, stimCellStr, colors, 1, gca, positiveSupport, 'vertical', 0);
  468. text(mean(xlim), 0.9 * getArrayElem(ylim, 2), ...
  469. sprintf('p = %1.1e', pVals(iAx + 4)), 'HorizontalAlignment', 'center');
  470. axis tight
  471. end
  472. % linkaxes(hSubAx(1:2), 'xy');
  473. % linkaxes(hSubAx(3:4), 'xy');
  474. ylabel(hSubAx(1), 'RF FWHM (deg)')
  475. hFig.Color = 'none';
  476. set(hSubAx, 'Color', 'none');
  477. arrayfun(@prettifyAxes, hSubAx)
  478. setFontForThesis(hSubAx, hFig);
  479. figFileName = ['CenterSurround-Width-Bean'];
  480. set(gcf,'renderer','painters')
  481. print(hFig, [figDir figFileName '.pdf'], '-dpdf')
  482. %% Export data for figure 7c-d
  483. dataCell = [offAmpCenter(:)', offAmpSurround(:)', onAmpCenter(:)', onAmpSurround(:)'];
  484. sheetName = ['Fig7cd'];
  485. startRow = 1;
  486. conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
  487. varNames = {[conditionStrCell{1} '_OFF_RF_Center_Amplitude_dF_F0'], ...
  488. [conditionStrCell{2} '_OFF_RF_Center_Amplitude_dF_F0'], ...
  489. [conditionStrCell{1} '_OFF_RF_Surround_Amplitude_dF_F0'], ...
  490. [conditionStrCell{2} '_OFF_RF_Surround_Amplitude_dF_F0'], ...
  491. [conditionStrCell{1} '_ON_RF_Center_Amplitude_dF_F0'], ...
  492. [conditionStrCell{2} '_ON_RF_Center_Amplitude_dF_F0'], ...
  493. [conditionStrCell{1} '_ON_RF_Surround_Amplitude_dF_F0'], ...
  494. [conditionStrCell{2} '_ON_RF_Surround_Amplitude_dF_F0']};
  495. capitalStr = 'ABCDEFGHIJKLMNOPQ';
  496. % Export response amplitudes
  497. for iVar = 1: numel(dataCell)
  498. DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
  499. writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
  500. 'Range', [capitalStr(iVar) num2str(startRow)]);
  501. end
  502. % Export FWHM
  503. dataCell = [offWidthCenter(:)', offWidthSurround(:)', onWidthCenter(:)', onWidthSurround(:)'];
  504. startRow = max(cellfun(@numel, offAmpCenter)) + 3;
  505. conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
  506. varNames = {[conditionStrCell{1} '_OFF_RF_Center_FWHM_deg'], ...
  507. [conditionStrCell{2} '_OFF_RF_Center_FWHM_deg'], ...
  508. [conditionStrCell{1} '_OFF_RF_Surround_FWHM_deg'], ...
  509. [conditionStrCell{2} '_OFF_RF_Surround_FWHM_deg'], ...
  510. [conditionStrCell{1} '_ON_RF_Center_FWHM_deg'], ...
  511. [conditionStrCell{2} '_ON_RF_Center_FWHM_deg'], ...
  512. [conditionStrCell{1} '_ON_RF_Surround_FWHM_deg'], ...
  513. [conditionStrCell{2} '_ON_RF_Surround_FWHM_deg']};
  514. capitalStr = 'ABCDEFGHIJKLMNOPQ';
  515. % Export response amplitudes
  516. for iVar = 1: numel(dataCell)
  517. DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
  518. writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
  519. 'Range', [capitalStr(iVar) num2str(startRow)]);
  520. end
  521. %% Export data for figure 7c-e stats
  522. dataCell = goodFlyIndsCell;
  523. sheetName = ['Fig7cd'];
  524. startRow = 2 * (max(cellfun(@numel, goodFlyIndsCell)) + 2) + 1;
  525. conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
  526. varNames = {[conditionStrCell{1} '_Fly_Index'], ...
  527. [conditionStrCell{2} '_Fly_Index']};
  528. capitalStr = 'ABCD';
  529. % Export response amplitudes
  530. for iVar = 1: numel(dataCell)
  531. DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
  532. writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
  533. 'Range', [capitalStr(iVar) num2str(startRow)]);
  534. end
  535. %% This is the single Gaussian fit version used for Supplementary figure 6.
  536. rSqThresh = -20;
  537. respIndThresh = 0.5;
  538. alignAll = true;
  539. [offX, offY, onX, onY] = deal(cell(size(StimTableNeuronStim([2, 1], 2: end), 1), 1));
  540. jTmp = 1;
  541. for iNeuron = [2, 1]
  542. StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
  543. [offX{jTmp}, offY{jTmp}, ...
  544. onX{jTmp}, onY{jTmp}] = getRFsAlignedWithFit(StimTableCell, ...
  545. rSqThresh, ...
  546. respIndThresh, alignAll);
  547. jTmp = jTmp + 1;
  548. end
  549. %% Inline from plotRFsAlignedWithFitGenotypes
  550. jTmp = 1;
  551. for iNeuron = [2, 1]
  552. StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
  553. [~, responseIndArrayCell, rSquareArrayCell, fitParamsCell, ...
  554. flyIndsCell, ~, ~] = getStimArraysForSameRois(StimTableCell);
  555. % Align to fit
  556. widthCell{jTmp} = cellfun(@(x) arrayfun(@(y) 2 * (2 * sqrt(log(2))) * y.c1, x, 'uni', 0), ...
  557. fitParamsCell, 'uni', 1); % The 2 factor is because of bar separation.
  558. validIndsCell{jTmp} = cellfun(@(x, y) x(:) > rSqThresh & y(:) > respIndThresh, ...
  559. rSquareArrayCell, responseIndArrayCell, 'uni', 0)';
  560. jTmp = jTmp + 1;
  561. end
  562. %%
  563. validOffInds = cellfun(@(x) x{1} & x{2}, validIndsCell, 'uni', 0);
  564. validOnInds = cellfun(@(x) x{3} & x{4}, validIndsCell, 'uni', 0);
  565. validAllInds = cellfun(@(x) x{1} & x{2} & x{3} & x{4}, validIndsCell, 'uni', 0);
  566. %%
  567. widthCell = cellfun(@(x) cat(1, x{:}), widthCell, 'uni', 0);
  568. %%
  569. xOFFWidth = cellfun(@(x, y) x(1, y), widthCell, validAllInds, 'uni', 0);
  570. yOFFWidth = cellfun(@(x, y) x(2, y), widthCell, validAllInds, 'uni', 0);
  571. xONWidth = cellfun(@(x, y) x(3, y), widthCell, validAllInds, 'uni', 0);
  572. yONWidth = cellfun(@(x, y) x(4, y), widthCell, validAllInds, 'uni', 0);
  573. allWidthsCell = {xOFFWidth, yOFFWidth, xONWidth, yONWidth};
  574. %%
  575. for iW = 1: numel(allWidthsCell)
  576. allWidthsCell{iW}(cellfun(@isempty, allWidthsCell{iW}, 'uni', 1)) = {pi};
  577. end
  578. %%
  579. offWidthCell = cellfun(@(x, y) (x + y) / 2, allWidthsCell{1}, allWidthsCell{2}, 'uni', 0);
  580. onWidthCell = cellfun(@(x, y) (x + y) / 2, allWidthsCell{2}, allWidthsCell{3}, 'uni', 0);
  581. %% Statistics
  582. % reviewer #1 analysis by fly
  583. % Number of ROIs differs across stimuli, but this has been accounted at
  584. % least for bars in getStimArraysForSameRois
  585. goodFlyIndsCell = cellfun(@(x, y) x(y), flyIndsCell(1: 2), validAllInds, 'uni', 0);
  586. nROIsPerFlyCell = cellfun(@(x) accumarray(x(:), 1), goodFlyIndsCell, 'uni', 0);
  587. % Delete flies with no rois to analyze.
  588. nROIsPerFlyCell{1}(nROIsPerFlyCell{1} == 0) = [];
  589. nROIsPerFlyCell{2}(nROIsPerFlyCell{2} == 0) = [];
  590. % generate continuos numbering.
  591. reorderedFlyIndsCell = cellfun(@(x) repelem(1: numel(x), x), nROIsPerFlyCell, 'uni', 0);
  592. % samples are cells, but we will permute all cells of a fly together.
  593. sample1 = 1: numel(nROIsPerFlyCell{1});
  594. sample2 = numel(sample1) + (1: numel(nROIsPerFlyCell{2}));
  595. %% Statistics for supplementary figure 6b
  596. paramsCell = {offWidthCell, onWidthCell};
  597. permutations = 10000;
  598. sidedness = 'both';
  599. plotresult = 1;
  600. pVals = nan(1, numel(paramsCell));
  601. obsDiffs = nan(1, numel(paramsCell));
  602. for iParam = 1: numel(paramsCell)
  603. % Permutation test needs row vectors
  604. paramsQuality = cellfun(@(x) x(:)', paramsCell{iParam}, 'uni', 0);
  605. [obsDiffs(iParam), pVals(iParam)] = nestetPermutationTest(sample1, sample2, ...
  606. paramsQuality, reorderedFlyIndsCell, ...
  607. permutations, sidedness, plotresult);
  608. end
  609. %% Plot supplementary figure 6b
  610. hFig = createPrintFig(15 * [1/4 1/3]);
  611. colors = brewermap(4, 'Set1');
  612. colors = colors([1, 2, 4, 3], :);
  613. hSubAx(1) = subplot(1, 2, 1);
  614. hWidthAx = beanPlot(offWidthCell, fliplr(cellTypeStr), colors, ...
  615. 1, gca, [0-eps 100], 'vertical', 0, 1, 0);
  616. text(mean(xlim), 0.9 * getArrayElem(ylim, 2), ...
  617. sprintf('p = %1.1e', pVals(1)), 'HorizontalAlignment', 'center');
  618. hSubAx(2) = subplot(1, 2, 2);
  619. hWidthAx = beanPlot(onWidthCell, fliplr(cellTypeStr), colors, ...
  620. 1, gca, [0-eps 100], 'vertical', 0, 1, 0);
  621. text(mean(xlim), 0.9 * getArrayElem(ylim, 2), ...
  622. sprintf('p = %1.1e', pVals(2)), 'HorizontalAlignment', 'center');
  623. axis(hSubAx(1), 'tight')
  624. axis(hSubAx(2), 'tight')
  625. % linkaxes(hSubAx, 'y');
  626. setFontForThesis(hWidthAx, gcf);
  627. arrayfun(@prettifyAxes, hWidthAx)
  628. % arrayfun(@offsetAxes, hWidthAx)
  629. figFileName = ['RF-FWHM-XYMean-bean'];
  630. print(hFig, [figDir figFileName '.pdf'], '-dpdf')
  631. %% Export data for supplementary figure 6b
  632. dataCell = [offWidthCell(:)', onWidthCell(:)'];
  633. sheetName = ['SuppFig6b'];
  634. startRow = 1;
  635. conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
  636. varNames = {[conditionStrCell{1} '_OFF_RF_Center_FWHM_deg'], ...
  637. [conditionStrCell{2} '_OFF_RF_Center_FWHM_deg'], ...
  638. [conditionStrCell{1} '_ON_RF_Center_FWHM_deg'], ...
  639. [conditionStrCell{2} '_ON_RF_Center_FWHM_deg']};
  640. capitalStr = 'ABCDEFGHIJKLMNOPQ';
  641. % Export response amplitudes
  642. for iVar = 1: numel(dataCell)
  643. DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
  644. writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
  645. 'Range', [capitalStr(iVar) num2str(startRow)]);
  646. end
  647. %% Min max orientation mean
  648. maxAmpCell = cellfun(@(x) max(x, [], 2), [offX offY onX onY], 'uni', 0);
  649. minAmpCell = cellfun(@(x) min(x, [], 2), [offX offY onX onY], 'uni', 0);
  650. offMaxCell = cellfun(@(x, y) (x + y) / 2, maxAmpCell(:, 1), maxAmpCell(:, 2), 'uni', 0);
  651. onMaxCell = cellfun(@(x, y) (x + y) / 2, maxAmpCell(:, 3), maxAmpCell(:, 4), 'uni', 0);
  652. offMinCell = cellfun(@(x, y) (x + y) / 2, minAmpCell(:, 1), minAmpCell(:, 2), 'uni', 0);
  653. onMinCell = cellfun(@(x, y) (x + y) / 2, minAmpCell(:, 3), minAmpCell(:, 4), 'uni', 0);
  654. %% Statistics for supplementary figure 6a
  655. paramsCell = {offMaxCell, offMinCell, onMinCell, onMaxCell};
  656. permutations = 10000;
  657. sidedness = 'both';
  658. plotresult = 1;
  659. pVals = nan(1, numel(paramsCell));
  660. obsDiffs = nan(1, numel(paramsCell));
  661. for iParam = 1: numel(paramsCell)
  662. % Permutation test needs row vectors
  663. paramsQuality = cellfun(@(x) x(:)', paramsCell{iParam}, 'uni', 0);
  664. [obsDiffs(iParam), pVals(iParam)] = nestetPermutationTest(sample1, sample2, ...
  665. paramsQuality, reorderedFlyIndsCell, ...
  666. permutations, sidedness, plotresult);
  667. end
  668. %% Plot supplementary figure 6a
  669. hFig = createPrintFig(15 * [1/2 1/3]);
  670. colors = getONOFFXYColors;
  671. colors = lines(2);
  672. stimCellStr = {'Control', 'CDM'};
  673. titleStrCell = {'OFF Max (center)', 'OFF Min (surround)', ...
  674. 'ON Min (center)', 'ON Max (surround)'};
  675. positiveSupport = [0-eps 1000];
  676. negativeSupport = [-1000 0-eps];
  677. supportCell = {negativeSupport, positiveSupport};
  678. supportInds = [2 1 1 2];
  679. for iParam = 1: numel(paramsCell)
  680. hSubAx(iParam) = subplot(1, numel(paramsCell), iParam);
  681. title(titleStrCell{iParam});
  682. beanPlot(paramsCell{iParam}, stimCellStr, colors, 1, gca, ...
  683. supportCell{supportInds(iParam)}, 'vertical', 0);
  684. text(mean(xlim), 0.9 * getArrayElem(ylim, 2), ...
  685. sprintf('p = %1.1e', pVals(iParam)), 'HorizontalAlignment', 'center');
  686. axis tight
  687. end
  688. % linkaxes(hSubAx(1:2), 'xy');
  689. % linkaxes(hSubAx(3:4), 'xy');
  690. ylabel(hSubAx(1), 'RF max (\DeltaF/F_0)')
  691. % hFig.Color = 'none';
  692. set(hSubAx, 'Color', 'none');
  693. arrayfun(@prettifyAxes, hSubAx)
  694. setFontForThesis(hSubAx, hFig);
  695. figFileName = ['MaxMinAmp-XYMean-Bean'];
  696. set(gcf,'renderer','painters')
  697. print(hFig, [figDir figFileName '.pdf'], '-dpdf')
  698. %% Export data for supplementary figure 6b
  699. dataCell = [offMaxCell(:)', offMinCell(:)', onMinCell(:)', onMaxCell(:)'];
  700. sheetName = ['SuppFig6a'];
  701. startRow = 1;
  702. conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
  703. varNames = {[conditionStrCell{1} '_OFF_RF_Maximum_Amplitude_dF_F0'], ...
  704. [conditionStrCell{2} '_OFF_RF_Maximum_Amplitude_dF_F0'], ...
  705. [conditionStrCell{1} '_OFF_RF_Minimum_Amplitude_dF_F0'], ...
  706. [conditionStrCell{2} '_OFF_RF_Minimum_Amplitude_dF_F0'], ...
  707. [conditionStrCell{1} '_ON_RF_Minimum_Amplitude_dF_F0'], ...
  708. [conditionStrCell{2} '_ON_RF_Minimum_Amplitude_dF_F0'], ...
  709. [conditionStrCell{1} '_ON_RF_Maximum_Amplitude_dF_F0'], ...
  710. [conditionStrCell{2} '_ON_RF_Maximum_Amplitude_dF_F0']};
  711. capitalStr = 'ABCDEFGHIJKLMNOPQ';
  712. % Export response amplitudes
  713. for iVar = 1: numel(dataCell)
  714. DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
  715. writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
  716. 'Range', [capitalStr(iVar) num2str(startRow)]);
  717. end
  718. %% Data for figure 7f-g
  719. for iNeuron = 1: size(StimTableNeuronStim, 1)
  720. allOnOffTracesCell = cat(1, StimTableNeuronStim{iNeuron, 1}.paddedEpochStacks{:});
  721. allOnOffTracesCell = cellfun(@(x) x', allOnOffTracesCell, 'uni', 0);
  722. allOffTraces = cat(1, allOnOffTracesCell{:, 1});
  723. allOnTraces = cat(1, allOnOffTracesCell{:, 2});
  724. allOnOffTracesArray = [allOffTraces allOnTraces];
  725. allNeuronsOnOffTraces{iNeuron} = allOnOffTracesArray;
  726. end
  727. %%
  728. % reviewer #1 analysis by fly
  729. flyIndsCdm = repelem(1:size(StimTableNeuronStim{1, 1}), ...
  730. cellfun(@(x) size(x, 2), StimTableNeuronStim{1, 1}.responseIndex));
  731. flyIndsNoCdm = repelem(1:size(StimTableNeuronStim{2, 1}), ...
  732. cellfun(@(x) size(x, 2), StimTableNeuronStim{2, 1}.responseIndex));
  733. %% Plot figure 7g
  734. cellTypeStr = {'CDM', 'noCDM'};
  735. desiredCellTypesStr = {'noCDM', 'CDM'};
  736. cellTypeInds = cellfun(@(x) find(strcmp(cellTypeStr, x)), desiredCellTypesStr);
  737. % cellTypeInds = 1:2;
  738. out = getONOFFParamsFromTable(StimTableNeuronStim);
  739. paramsCell = {out.offMean, out.onMean, out.responseInd};
  740. supportCell = {'unbounded', 'unbounded', [-1.01 1.01], [-1 1], [0 2], [-eps 2.01], [0 1]};
  741. paramsNames = {'Mean OFF', 'Mean ON', 'Polarity index', 'Sustenance index', ...
  742. 'Half-rise time (s)', 'Time to extreme (s)', 'Response quality'};
  743. qualityThreshold = 0.5;
  744. paramsCell = cellfun(@(x) x(cellTypeInds), paramsCell, 'uni', 0);
  745. colors = repmat([1 0.75 1] * 0.25, numel(cellTypeInds), 1);
  746. colors = lines(2);
  747. hFig = createPrintFig(15 * [1/3 1]);
  748. for iParam = 1: numel(paramsCell) - 1
  749. hSubAx(iParam) = subplot(2, 1, iParam);
  750. title(paramsNames{iParam});
  751. paramsQuality = cellfun(@(x, y) x(y > qualityThreshold), paramsCell{iParam}, paramsCell{end}, 'uni', 0);
  752. paramsToPlot = paramsQuality;
  753. % paramsToPlot = paramsCell{iParam};
  754. if 1
  755. beanPlot(paramsToPlot, desiredCellTypesStr, colors, 1, gca, ...
  756. supportCell{iParam}, 'vertical', 0, 0, 1);
  757. arrayfun(@prettifyAxes, gca);
  758. arrayfun(@offsetAxes, gca);
  759. setFontForThesis(gca, hFig)
  760. else
  761. boxInd = repelem(1:numel(paramsToPlot), cellfun(@numel, paramsToPlot));
  762. boxplot(cat(2, paramsToPlot{:}), boxInd, 'Notch', 'on')
  763. end
  764. if numel(paramsToPlot) > 1
  765. for jNeuron = 2: numel(paramsToPlot)
  766. [pVal(iParam, jNeuron - 1), ...
  767. obsDiff(iParam, jNeuron - 1)] = permutationTest(paramsToPlot{1}, ...
  768. paramsToPlot{jNeuron}, ...
  769. 50000, 'plotresult', 0);
  770. end
  771. end
  772. end
  773. setFontForThesis(gca, hFig)
  774. figFileName = ['OnOff-Params-'];
  775. set(gcf,'renderer','painters')
  776. % set(gcf, 'PaperPositionMode', 'auto');
  777. print(hFig, [figDir figFileName '.pdf'], '-dpdf')
  778. print(hFig, [figDir figFileName '.png'], '-dpng', '-r300')
  779. %% Export data for figure 7g
  780. paramsQuality = arrayfun(@(z) cellfun(@(x, y) x(y > qualityThreshold), paramsCell{z}, paramsCell{end}, 'uni', 0), 1:2,'uni', 0);
  781. dataCell = [paramsQuality{1}(:)', paramsQuality{2}(:)'];
  782. sheetName = ['Fig7g'];
  783. startRow = 1;
  784. conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
  785. varNames = {[conditionStrCell{1} '_OFF_Fullfield_Mean_Amplitude_dF_F0'], ...
  786. [conditionStrCell{2} '_OFF_Fullfield_Mean_Amplitude_dF_F0'], ...
  787. [conditionStrCell{1} '_ON_Fullfield_Mean_Amplitude_dF_F0'], ...
  788. [conditionStrCell{2} '_ON_Fullfield_Mean_Amplitude_dF_F0']};
  789. capitalStr = 'ABCDEFGHIJKLMNOPQ';
  790. % Export response amplitudes
  791. for iVar = 1: numel(dataCell)
  792. DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
  793. writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
  794. 'Range', [capitalStr(iVar) num2str(startRow)]);
  795. end
  796. %% plot figure 7f
  797. onOffTraces = cellfun(@(x, y) x(y > qualityThreshold, :), ...
  798. allNeuronsOnOffTraces(cellTypeInds), out.responseInd(cellTypeInds), 'uni', 0);
  799. colors = repmat([1 0.75 1] * 0.25, numel(cellTypeStr), 1);
  800. % colors = brewermap(numel(desiredCellTypesStr), 'Dark2');
  801. colors = brewermap([], 'Dark2');
  802. colors = colors([8, (1: numel(cellTypeStr) - 1) + 3*0], :);
  803. % For cdm
  804. colors = lines(2);
  805. plotONOFFTracesFromCell(onOffTraces, desiredCellTypesStr, figDir, colors)
  806. %% Export data for figure 7f
  807. dataFileName = [figDir filesep 'Fig7.xls'];
  808. timePoints = (1: numel(onOffTraces{1}(1, :))) / 10; % Data interpolated at 10Hz
  809. sheetName = 'Fig7f';
  810. stimNameStrCell = {'Tm9_Control_OFFON_fullfield_time_seconds', ...
  811. 'Tm9_CDM_OFFON_fullfield_time_seconds'};
  812. startRow = 1;
  813. for iStim = 1: 2
  814. DataTable = createRFsTable(onOffTraces{iStim}, stimNameStrCell{iStim});
  815. DataTable(:, 1).(stimNameStrCell{iStim}) = timePoints(:);
  816. writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
  817. startRow = startRow + size(DataTable, 1) + 2;
  818. end
  819. %% More sophisticated test for reviewer #1
  820. flyIndsCell = {flyIndsNoCdm flyIndsCdm};
  821. goodFlyIndsCell = cellfun(@(x, y) x(y > qualityThreshold), flyIndsCell, paramsCell{end}, 'uni', 0);
  822. flyGroups = [goodFlyIndsCell{1} max(goodFlyIndsCell{1}) + goodFlyIndsCell{2}];
  823. genotypeGroups = [repmat({'control'}, 1, numel(goodFlyIndsCell{1})) ...
  824. repmat({'cdm'}, 1, numel(goodFlyIndsCell{2}))];
  825. groups = {genotypeGroups, flyGroups};
  826. % I think the fly level is a random effect, right?
  827. pValsCell = cell(1, numel(paramsCell));
  828. for iParam = 1: numel(paramsCell)
  829. paramsQuality = cellfun(@(x, y) x(y > qualityThreshold), paramsCell{iParam}, paramsCell{end}, 'uni', 0);
  830. [pValsCell{iParam}, ~, statsCell{iParam}] = ...
  831. anovan(cat(2, paramsQuality{:}), ...
  832. groups, ...
  833. 'nested',[0 0; 1 0], ...
  834. 'varnames', {'Genotype', 'Fly'}, ...
  835. 'random', [2], 'display', 'on');
  836. end
  837. %%
  838. hFig = createPrintFig(10 * [1 1/2]);
  839. hIm = imagesc([], [1, 2], cat(2, pValsCell{:}));
  840. colormap([plasma; 0.5 0.5 0.5]);
  841. caxis([0 0.051]);
  842. set(gca, 'XTick', 1: numel(paramsNames), 'XTickLabels', paramsNames, 'XTickLabelRotation', 90);
  843. set(gca, 'YTick', [1 2], 'YTickLabels', {'Pooling across cells', 'Pooling across flies'}, ...
  844. 'XTickLabelRotation', 30)
  845. title('p-values from nested anova FFF');
  846. colorbar;
  847. arrayfun(@prettifyAxes, gca)
  848. setFontForThesis(gca, hFig);
  849. figFileName = ['NestedAnova-pValues-ONOFF-FFF-'];
  850. set(gcf,'renderer','painters')
  851. print(hFig, [figDir figFileName '.pdf'], '-dpdf');
  852. %% Still do the nested permutation
  853. nROIsPerFlyCell = cellfun(@(x) accumarray(x(:), 1), goodFlyIndsCell, 'uni', 0);
  854. % Delete flies with no rois to analyze.
  855. nROIsPerFlyCell{1}(nROIsPerFlyCell{1} == 0) = [];
  856. nROIsPerFlyCell{2}(nROIsPerFlyCell{2} == 0) = [];
  857. % generate continuos numbering.
  858. reorderedFlyIndsCell = cellfun(@(x) repelem(1: numel(x), x), nROIsPerFlyCell, 'uni', 0);
  859. % samples are cells, but we will permute all cells of a fly together.
  860. sample1 = 1: numel(nROIsPerFlyCell{1});
  861. sample2 = numel(sample1) + (1: numel(nROIsPerFlyCell{2}));
  862. %% Statistics for figure 7g
  863. permutations = 10000;
  864. sidedness = 'both';
  865. plotresult = 1;
  866. pVals = nan(1, numel(paramsCell));
  867. obsDiffs = nan(1, numel(paramsCell));
  868. for iParam = 1: numel(paramsCell)
  869. paramsQuality = cellfun(@(x, y) x(y > qualityThreshold), paramsCell{iParam}, paramsCell{end}, 'uni', 0);
  870. [obsDiffs(iParam), pVals(iParam)] = nestetPermutationTest(sample1, sample2, ...
  871. paramsQuality, reorderedFlyIndsCell, ...
  872. permutations, sidedness, plotresult);
  873. end
  874. %% Plot the pvalues as bars
  875. hFig = createPrintFig(15 * [1 3/4]);
  876. bar(pVals, 'FaceColor', [1 1 1] * 0.15, 'EdgeColor', 'none');
  877. arrayfun(@(x) text(x, 0.02 + pVals(x), num2str(pVals(x), '%1.2e'), ...
  878. 'HorizontalAlignment', 'center'), 1: numel(pVals))
  879. ylabel('p-value')
  880. hRefLine = refline(0, 0.05);
  881. hRefLine.Color = [1 1 1] * 0.5;
  882. hRefLine.LineStyle = ':';
  883. set(gca, 'XTick', 1: numel(paramsNames), ...
  884. 'XTickLabels', paramsNames, 'XTickLabelRotation', 30);
  885. set(gca, 'YTick', [0, 0.01, 0.05, 0.1, 0.2]);
  886. title('p-values from nested permutation test');
  887. arrayfun(@prettifyAxes, gca)
  888. setFontForThesis(gca, hFig);
  889. figFileName = ['ONOFF-NestedPermutations-pValues'];
  890. set(gcf,'renderer','painters')
  891. print(hFig, [figDir figFileName '.pdf'], '-dpdf');
  892. %% Export data for figure 7g stats
  893. dataCell = goodFlyIndsCell;
  894. sheetName = ['Fig7g'];
  895. startRow = 1;
  896. conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
  897. varNames = {[conditionStrCell{1} '_Fly_Index'], ...
  898. [conditionStrCell{2} '_Fly_Index']};
  899. capitalStr = 'EFGHIJKLMNOPQ';
  900. % Export response amplitudes
  901. for iVar = 1: numel(dataCell)
  902. DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
  903. writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
  904. 'Range', [capitalStr(iVar) num2str(startRow)]);
  905. end
  906. %%
  907. function [Amp, Pos, Sigma, RSq, varargout] = parseFitParams(fitStructCellArray)
  908. TmpGofCell = cellfun(@(y) arrayfun(@(x) x.gof.rsquare, y(1: end-1)), fitStructCellArray, 'uni', 0);
  909. if isstruct(fitStructCellArray{1, 1}(1).fit)
  910. fitParams = fieldnames(fitStructCellArray{1, 1}(1).fit);
  911. ParamCell = cell(1, numel(fitParams));
  912. for iField = 1: numel(fitParams)
  913. % Delete the last ROI from background
  914. fieldCell = cellfun(@(y) arrayfun(@(x) x.fit.(fitParams{iField}), y(1: end-1)), fitStructCellArray, 'uni', 0);
  915. ParamCell{iField} = cell(1, size(fieldCell, 2));
  916. for jFly = 1: size(ParamCell{iField}, 2)
  917. ParamCell{iField}{jFly} = cat(1, fieldCell{:, jFly});
  918. end
  919. end
  920. else
  921. % Params are given as an array.
  922. ParamCell = cell(1, numel(fitStructCellArray{1, 1}(1).fit));
  923. for iField = 1: numel(fitStructCellArray{1, 1}(1).fit)
  924. % Delete the last ROI from background
  925. fieldCell = cellfun(@(y) arrayfun(@(x) x.fit(iField), y(1: end-1))', fitStructCellArray, 'uni', 0);
  926. fieldCell = squeeze(fieldCell);
  927. ParamCell{iField} = cell(1, size(fieldCell, 2));
  928. for jFly = 1: size(ParamCell{iField}, 2)
  929. ParamCell{iField}{jFly} = cat(1, fieldCell{:, jFly});
  930. end
  931. end
  932. TmpGofCell = cellfun(@(x) x', TmpGofCell, 'uni', 0);
  933. end
  934. for jFly = 1: size(TmpGofCell, 2)
  935. RSq{jFly} = cat(1, TmpGofCell{:, jFly});
  936. end
  937. if numel(ParamCell) > 5
  938. [Amp.center, Pos, Sigma.center, Amp.surround, Sigma.surround, baseline] = deal(ParamCell{:});
  939. else
  940. [Amp.center, Pos, Sigma.center, Amp.surround, Sigma.surround] = deal(ParamCell{:});
  941. end
  942. if nargout > 4
  943. if numel(ParamCell) > 5
  944. varargout{1} = baseline;
  945. else
  946. varargout{1} = nan;
  947. end
  948. end
  949. end