diff --git a/+helper/bonf_holm.m b/+helper/bonf_holm.m new file mode 100644 index 0000000..e50b0d3 --- /dev/null +++ b/+helper/bonf_holm.m @@ -0,0 +1,113 @@ + +% Bonferroni-Holm (1979) correction for multiple comparisons. This is a +% sequentially rejective version of the simple Bonferroni correction for multiple +% comparisons and strongly controls the family-wise error rate at level alpha. +% +% It works as follows: +% 1) All p-values are sorted in order of smallest to largest. m is the +% number p-values. +% 2) If the 1st p-value is greater than or equal to alpha/m, the procedure +% is stopped and no p-values are significant. Otherwise, go on. +% 3) The 1st p-value is declared significant and now the second p-value is +% compared to alpha/(m-1). If the 2nd p-value is greater than or equal +% to alpha/(m-1), the procedure is stopped and no further p-values are +% significant. Otherwise, go on. +% 4) Et cetera. +% +% As stated by Holm (1979) "Except in trivial non-interesting cases the +% sequentially rejective Bonferroni test has strictly larger probability of +% rejecting false hypotheses and thus it ought to replace the classical +% Bonferroni test at all instants where the latter usually is applied." +% +% +% function [corrected_p, h]=bonf_holm(pvalues,alpha) +% +% Required Inputs: +% pvalues - A vector or matrix of p-values. If pvalues is a matrix, it can +% be of any dimensionality (e.g., 2D, 3D, etc...). +% +% Optional Input: +% alpha - The desired family-wise alpha level (i.e., the probability of +% rejecting one of more null hypotheses when all null hypotheses are +% really true). {default: 0.05} +% +% Output: +% corrected_p - Bonferroni-Holm adjusted p-values. Any adjusted p-values +% less than alpha are significant (i.e., that null hypothesis +% is rejected). The adjusted value of the smallest p-value +% is p*m. The ith smallest adjusted p-value is the max of +% p(i)*(m-i+1) or adj_p(i-1). Note, corrected p-values can +% be greater than 1. +% h - A binary vector or matrix of the same dimensionality as +% pvalues. If the ith element of h is 1, then the ith p-value +% of pvalues is significant. If the ith element of h is 0, then +% the ith p-value of pvalues is NOT significant. +% +% Example: +% >>p=[.56 .22 .001 .04 .01]; %five hypothetical p-values +% >>[cor_p, h]=bonf_holm(p,.05) +% cor_p = +% 0.5600 0.4400 0.0050 0.1200 0.0400 +% h = +% 0 0 1 0 1 +% +% Conclusion: the third and fifth p-values are significant, but not the +% remaining three. +% +% Reference: +% Holm, S. (1979) A simple sequentially rejective multiple test procedure. +% Scandinavian Journal of Statistics. 6, 65-70. +% +% +% For a review on contemporary techniques for correcting for multiple +% comparisons that are often more powerful than Bonferroni-Holm see: +% +% Groppe, D.M., Urbach, T.P., & Kutas, M. (2011) Mass univariate analysis +% of event-related brain potentials/fields I: A critical tutorial review. +% Psychophysiology, 48(12) pp. 1711-1725, DOI: 10.1111/j.1469-8986.2011.01273.x +% http://www.cogsci.ucsd.edu/~dgroppe/PUBLICATIONS/mass_uni_preprint1.pdf +% +% +% Author: +% David M. Groppe +% Kutaslab +% Dept. of Cognitive Science +% University of California, San Diego +% March 24, 2010 +function [corrected_p, h]=bonf_holm(pvalues,alpha) +if nargin<1 + error('You need to provide a vector or matrix of p-values.'); +else + if ~isempty(find(pvalues<0,1)), + error('Some p-values are less than 0.'); + elseif ~isempty(find(pvalues>1,1)), + fprintf('WARNING: Some uncorrected p-values are greater than 1.\n'); + end +end +if nargin<2 + alpha=.05; +elseif alpha<=0 + error('Alpha must be greater than 0.'); +elseif alpha>=1 + error('Alpha must be less than 1.'); +end +s=size(pvalues); +if isvector(pvalues) + if size(pvalues,1)>1 + pvalues=pvalues'; + end + [sorted_p sort_ids]=sort(pvalues); +else + [sorted_p sort_ids]=sort(reshape(pvalues,1,prod(s))); +end +[dummy, unsort_ids]=sort(sort_ids); %indices to return sorted_p to pvalues order +sorted_p(isnan(sorted_p)) = []; +m=length(sorted_p); %number of tests +mult_fac=m:-1:1; +cor_p_sorted=sorted_p.*mult_fac; +cor_p_sorted(2:m)=max([cor_p_sorted(1:m-1); cor_p_sorted(2:m)]); %Bonferroni-Holm adjusted p-value +cor_p_sorted = [cor_p_sorted NaN(1,length(sort_ids)-length(cor_p_sorted))]; +corrected_p=cor_p_sorted(unsort_ids); +corrected_p=reshape(corrected_p,s); +h=corrected_p lickFreqThr)); + +lickErrorsInd = union(shortLickInd,tooFastLickInd); +lickStartInd = setdiff(lickStartInd,lickErrorsInd); + +lickEvents(1:numel(lickSignal)) = zeros; +lickEvents(lickStartInd) = 1; + +lickFrequency = helper.ifreq(lickEvents,samplingRate); +lickFrequencyDs = lickFrequency(sData.daqdata.frameIndex(1):323:end); + + +% lickFrequencyDs =[]; +end \ No newline at end of file diff --git a/+helper/getStats.m b/+helper/getStats.m new file mode 100644 index 0000000..4b194da --- /dev/null +++ b/+helper/getStats.m @@ -0,0 +1,58 @@ +function roiStat = getStats(sData, CH) +%getRoiActivityStats Return struct of different roi activity statistics +% +% roiStat = getRoiActivityStats(sData) Returns struct with statistics +% about the roi signal for all rois in sData. +% +% roiStat = getRoiActivityStats(sData, CH) calculates the roiStat on the +% specified channel number CH. CH is an integer. The default value is 2. +% +% Output, roiStat contains the following fields: +% peakDff : Peak dff of the roi time series. +% signalToNoise : Signal to noise +% activityLevel : Fraction of time where activity is above noise +% + +% % pstr = getFilePath(sData.sessionID, 'roiStat'); +% % if exist(pstr, 'file') +% % roiStat = loaddata(sData.sessionID, 'roiStat'); +% % return +% % end + + if nargin < 2 + CH = 2; % Most of us are imaging on ch2? + end + + dff = double(squeeze(sData.imdata.roiSignals(CH).dff)); + [nRois, nSamples] = size(dff); + + + % Get max DFF of all rois. + peakDff = max(dff, [], 2); + + + % Get SNR of all Rois. + signalToNoise = zeros(nRois, 1); + noiseLevel = zeros(nRois, 1); + for i = 1:nRois + noiseLevel(i) = real(GetSn(dff(i, :))); + signalToNoise(i) = snr(dff(i, :), ones(1,nSamples) * noiseLevel(i)); + end + + signalToNoise = helper.calculateSNR(sData); + + % Get fraction of time above noise level + %dffSmooth = okada(dff, 2); % Smooth with okada filter + dffSmooth = smoothdata(dff, 2, 'movmean', 7); % Smooth again with movmean + + isHigh = dffSmooth > noiseLevel; + + activityLevel = sum(isHigh, 2) ./ nSamples; + + roiStat = struct('peakDff', peakDff, ... + 'signalToNoise', signalToNoise, .... + 'activityLevel', activityLevel); + +% savedata(sData.sessionID, struct('roiStat', roiStat)) + +end \ No newline at end of file diff --git a/+helper/ifreq.m b/+helper/ifreq.m new file mode 100644 index 0000000..b63ce0e --- /dev/null +++ b/+helper/ifreq.m @@ -0,0 +1,16 @@ +function iFrequency = ifreq(signal,fs) + +% returns the instantaneous frequency of a discrete signal (containing 0 and 1) with fs sampling rate + +iFrequency = nan(size(signal)); +if size(signal,1) > size(signal,2) + signal = signal'; +end +events = find(signal == 1); +freques = fs./diff(events); +freques = [freques 0]; +iFrequency(signal == 1) = freques; +iFrequency(1) = 0; +iFrequency = fillmissing(iFrequency,'previous'); + +end \ No newline at end of file diff --git a/+helper/splitInTrialsPos.m b/+helper/splitInTrialsPos.m new file mode 100644 index 0000000..8455488 --- /dev/null +++ b/+helper/splitInTrialsPos.m @@ -0,0 +1,13 @@ +function binnedTrialedSingleROI = splitInTrialsPos(signal, sData,idx) + +for i = 1:numel(sData.trials.trialStart)-1 + isInTrial = sData.trials.trialLength == i; + trialSignal = signal(idx,isInTrial); + + for jj = 1:78 + idxInTrialBin = sData.behavior.wheelPosDsBinned(isInTrial) == jj; + binnedTrialedSingleROI(i,jj) = nanmean(trialSignal(idxInTrialBin)); + end +end + +end \ No newline at end of file diff --git a/+plot/createHeatMap.m b/+plot/createHeatMap.m new file mode 100644 index 0000000..6c6989f --- /dev/null +++ b/+plot/createHeatMap.m @@ -0,0 +1,23 @@ + +function matrix = createHeatMap(lap, pos, behavior) + + % make all the same length + minLength = min([length(lap),length(pos), length(behavior)]); + lap = lap(1:minLength); + pos = pos(1:minLength); + behavior = behavior(1:minLength); + + heatmap = zeros(max(lap)+2,ceil(max(pos))); + occupancymap = zeros(max(lap)+2,ceil(max(pos))); + offset = min(lap); + for x = 1:length(lap)%pos) + x_pos = round(pos(x))+1; + if x_pos > size(heatmap,2); x_pos = size(heatmap,2); end + y_pos = lap(x)+1+(-offset); + heatmap(y_pos,x_pos) = heatmap(y_pos,x_pos) + behavior(x); + occupancymap(y_pos,x_pos) = occupancymap(y_pos,x_pos) + 1; + end + + matrix = heatmap./occupancymap; + +end diff --git a/+plot/decodingVsNoCells.m b/+plot/decodingVsNoCells.m new file mode 100644 index 0000000..0f71988 --- /dev/null +++ b/+plot/decodingVsNoCells.m @@ -0,0 +1,55 @@ +rmse_all = NaN(length(data(area).sessionIDs),1); +cell_no = cell(length(data(area).sessionIDs),1); +cellsSpaced = NaN(length(data(area).sessionIDs),1); +lengthCells = NaN(length(data(area).sessionIDs),1); + +f =1; +for g = 1: length(data(area).sessionIDs) + + % read out decoding error + rmse_session = NaN(1,length(mData(area,f).d_data)); + if length(mData(area,g).d_data)>1 + for i = 1:length(mData(area,g).d_data) + cell_no{f}(i) = length(mData(area,g).d_data(i).iter{1, mData(area,g).d_data(i).medianTestFold}.cellIdx); + rmse_session(i) = mData(area,g).d_data(i).meanRMSE_test; + end + + % prepare a matrix that contains all decoding errors + maxLength = max([length(rmse_session),size(rmse_all,2)]); + rmse_all = [rmse_all NaN(size(rmse_all,1), maxLength-size(rmse_all,2))]; + rmse_session = [rmse_session NaN(1,maxLength-length(rmse_session))]; + + rmse_all(f,:) = rmse_session; + + cellsSpaced(f) = 10*round(cell_no{f}(2)/10); + lengthCells(f) = length(cell_no{f}); + f =f+1; + end + +end + + +% make sure that the spacing of no of cells used for decoding is the +% same for all sessions (starts at 10 and is then spaced with 20) +% this is necessary because for some sessions we calculated the decoding erroor with spacing of 10 +% cells +if ~isempty(find(cellsSpaced == 30, 1)) + cellsSpacedShort = find(cellsSpaced == 20); + lengthCells(cellsSpacedShort) = ceil(lengthCells(cellsSpacedShort)/2); + for f = 1:length(cellsSpacedShort) + rmse_all(cellsSpacedShort(f),2:2:end) = NaN; + end +end + +[maxCells, maxCellIdx]= max(lengthCells); +cellNoSum = cell_no{maxCellIdx}; + +% find idx from which there is only one session left, these will not be use +noNan = NaN(size(rmse_all,2),1); +for m = 1:size(rmse_all,2) + noNan(m) = length(find(isnan(rmse_all(:,m)))); +end + +idx_singleSesssion = find(noNan == size(rmse_all,1)-1); +rmse_all(:,idx_singleSesssion) = []; +cellNoSum(idx_singleSesssion) = []; diff --git a/+plot/fraction_tuning.m b/+plot/fraction_tuning.m new file mode 100644 index 0000000..74e2149 --- /dev/null +++ b/+plot/fraction_tuning.m @@ -0,0 +1,90 @@ + +col = [215,223,225; + 175 200 211; + 80 106 139; + 142 175 174; + 230 151 64; + 218 57 38]./255; + + +prct_all = cell(1); +prct_all_threeCategories = cell(1); +prct_Motor = cell(1); + +for l = 1:length(areas) + + for f = 1: length(data(areas(l)).sessionIDs) + + cue_cells = 0; + place_cells = 0; + speed_cells = 0; + acceleration_cells = 0; + lick_cells = 0; + not_tuned_cells = 0; + + for i = 1:length(mData(areas(l),f).LNP_results) + + if ~isempty(find(mData(areas(l),f).LNP_results(i).selected_model == 5, 1)) + cue_cells = cue_cells+1; + else + if ~isempty(find(mData(areas(l),f).LNP_results(i).selected_model == 1, 1)) + place_cells = place_cells+1; end + if ~isempty(find(mData(areas(l),f).LNP_results(i).selected_model == 2, 1)) + speed_cells =speed_cells+1; end + if ~isempty(find(mData(areas(l),f).LNP_results(i).selected_model == 3, 1)) + acceleration_cells = acceleration_cells+1; end + if ~isempty(find(mData(areas(l),f).LNP_results(i).selected_model == 4, 1)) + lick_cells= lick_cells+1; end + if isnan(mData(areas(l),f).LNP_results(i).selected_model) + not_tuned_cells = not_tuned_cells+1; end + end + end + + + sumval = not_tuned_cells+place_cells+speed_cells+acceleration_cells+lick_cells+cue_cells; + prct_session = 100*[not_tuned_cells place_cells speed_cells acceleration_cells lick_cells cue_cells]./length(mData(areas(l),f).LNP_results); + prct_all{l}(f,:) = prct_session; + + + % we are mainly interested in no of place tuned cells: therefor we + % merge all other categories + prct_session_threeCategories = 100*[not_tuned_cells place_cells speed_cells+acceleration_cells+lick_cells cue_cells]./length(mData(areas(l),f).LNP_results); + prct_all_threeCategories{l}(f,:) = prct_session_threeCategories; + prct_Motor{l}(f,:) = 100*[speed_cells acceleration_cells lick_cells]./length(mData(areas(l),f).LNP_results); + + end + + fig = figure(); + subplot(1,2,1) + meanvals = nanmean(prct_all_threeCategories{l},1); + meanstd = nanstd(prct_all_threeCategories{l},1); + for m = 2: length(meanvals) + hold on + scatter(m-1,meanvals(m), 70,col(m,:)); + errorbar(m-1,meanvals(m),meanstd(m)/sqrt(length(data(l).sessionIDs)),'Color',col(m,:),'LineWidth',1) + end + + xlim([0.5 3.5]) + xticks([1:3]) + xticklabels({'P','M','C'}) + ylim([0 50]) + ylabel('Tuning (%)') + set(gca,'FontName','Arial','FontSize',12) + title(data(areas(l)).area) + + subplot(1,2,2) + meanvals = nanmean(prct_Motor{l},1); + meanstd = nanstd(prct_Motor{l},1); + for m = 1: length(meanvals) + hold on + scatter(m,meanvals(m), 70,col(3+m,:)); + errorbar(m,meanvals(m),meanstd(m)/sqrt(length(data(l).sessionIDs)),'Color',col(3+m,:),'LineWidth',1) + end + ylim([0 15]) + xlim([0.5 3.5]) + xticks([1:3]) + xticklabels({'S','A','L'}) + ylabel('Tuning (%)') + set(gca,'FontName','Arial','FontSize',12) + title(data(areas(l)).area) +end \ No newline at end of file diff --git a/+plot/peakAlignedSessionMean.m b/+plot/peakAlignedSessionMean.m new file mode 100644 index 0000000..6344dd3 --- /dev/null +++ b/+plot/peakAlignedSessionMean.m @@ -0,0 +1,57 @@ +function peakAlignedSessionMean(meanSignal_area) + + +col = [122 172 210 + 80 106 139 + 217 85 88 + 129 89 162 + 180 151 94 + 179 205 142 + 0 0 0 + 0 0 0]./255; + + fig = figure(); + for i = 1:7 + subplot(3,3,i) + % for m = 1:length(file{i}) + normedSignalAll = (meanSignal_area{i} - min(meanSignal_area{i},[],2))./(max(meanSignal_area{i},[],2)-min(meanSignal_area{i},[],2)); + + meanSequence = nanmean(meanSignal_area{i},1); + plot(meanSequence,'Color',col(i,:),'LineWidth',1) + yval = meanSignal_area{i}; + xval = 1:157; + patch(gca,[xval fliplr(xval)],[nanmean(yval)+nanstd(yval)./sqrt(size(yval,1)) fliplr(nanmean(yval)-nanstd(yval)./sqrt(size(yval,1)))],col(i,:),'linestyle','none','FaceAlpha', 0.2); + + hold on + % end + +% ax = gca; +% ax.FontSize = 12; +% ax.FontName = 'Arial'; +% ax.XColor = [0 0 0]; +% ax.YColor = [0 0 0]; +% ax.XLabel.Color = [0 0 0]; +% ax.YLabel.Color = [0 0 0]; +% xticks([1 20 40 60 78]) +% xticklabels([0 40 60 80 157]) +% xlim([1 78]) +% xlabel('Position (cm)') +% ylabel('Deconv. DFF') +% box off +% title(area{i}) + end + + +ax = gca; +ax.FontSize = 12; +ax.FontName = 'Arial'; +ax.XColor = [0 0 0]; +ax.YColor = [0 0 0]; +ax.XLabel.Color = [0 0 0]; +ax.YLabel.Color = [0 0 0]; +xticks([1 20 40 60 78]) +xticklabels([0 40 60 80 157]) +xlim([1 78]) +xlabel('Position (cm)') +ylabel('Deconv. DFF') +box off \ No newline at end of file diff --git a/+plot/plotRealVsDecodedPosition.m b/+plot/plotRealVsDecodedPosition.m new file mode 100644 index 0000000..c405398 --- /dev/null +++ b/+plot/plotRealVsDecodedPosition.m @@ -0,0 +1,22 @@ +function plotRealVsDecodedPosition(mData,area,session) + col = [122 172 210 + 80 106 139 + 217 85 88 + 129 89 162 + 180 151 94 + 179 205 142]./255; + + real_position = mData(area ,session).d_data(length(mData(area ,session).d_data)).iter{1,mData(area ,session).d_data(length(mData(area ,session).d_data)).medianTestFold}.realPos_test; + predicted_position = mData(area ,session).d_data(length(mData(area ,session).d_data)).iter{1,mData(area ,session).d_data(length(mData(area ,session).d_data)).medianTestFold}.predPos_test; + + figure() + plot(real_position, 'LineWidth', 1.5, 'Color', 'k'); hold on + plot(predicted_position, 'LineWidth',1.5, 'Color', col(area,:)); + xlim( [100 700]) + xticks([]); + yticks([0 156]) + yticklabels({'0','157'}); + box 'off' + set(gca, 'FontSize',12, 'FontName', 'Arial') + axis('off') +end \ No newline at end of file diff --git a/+plot/plotSortedSequence.m b/+plot/plotSortedSequence.m new file mode 100644 index 0000000..739770e --- /dev/null +++ b/+plot/plotSortedSequence.m @@ -0,0 +1,49 @@ + + +for k =1:length(areas) + meanSignal_area(k).deconv = []; + for j = 1:length(data(areas(k)).sessionIDs) + for f = 1:length(mData(areas(k),j).LNP_results) + if ~isempty(find(mData(areas(k),j).LNP_results(f).selected_model == 1)) && isempty(find(mData(areas(k),j).LNP_results(f).selected_model == 5)) + meanSignal_area(k).deconv = [ meanSignal_area(k).deconv; rmaps(areas(k),j).deconv(f,:)]; + end + end + end + + normedSignal = (meanSignal_area(k).deconv - min(meanSignal_area(k).deconv, [], 2))./(max(meanSignal_area(k).deconv, [], 2) - min(meanSignal_area(k).deconv, [], 2)); + + [~, idx] = max(normedSignal, [], 2); % Returns 2 arrays. Val is the max value of each ROI, idx is the column index of the max value of each ROI. + idxMatrix = horzcat(idx, normedSignal); % Appends idx column array to the ROI matrix. + sortedMatrix = sortrows(idxMatrix, 1); % Sorts ROI matrix by the index of their max values in ascending order. + sortedMatrix(: , 1) = []; % Removes indexes from first column of sortedMatrix. + + fig = figure(); + imagesc(sortedMatrix) + c = colorbar; + colormap(jet); + c.Label.String = 'Deconv Signal'; + c.Label.FontSize = 11; + c.Label.FontUnits = 'normalized'; + c.TickDirection = 'out'; + caxis([0 1]); % Sets limit for color map, left is min (black) and right is max (white). + xlabel('Position (cm)'); + title(data(areas(k)).area) + xticks(gca,[1 25 50 78]) + xticklabels(gca,{'0' ,'50','100','157'}) + xlim([1 79]) + set(gca, 'FontSize',18, 'FontName', 'Arial') + + hold on + cue1 = [42 47 122 128]./2; + cue2 = [63 68 100 106]./2; + limsY = get(gca,'YLim'); + + for i = 1: length(cue1) + plot(gca, [cue1(i) cue1(i)],limsY, 'LineStyle','--','Color', [1 1 1], 'LineWidth', 2); + end + + for i = 1: length(cue2) + plot(gca, [cue2(i) cue2(i)],limsY, 'LineStyle','--','Color', [1 1 1], 'LineWidth', 2); + end + +end \ No newline at end of file diff --git a/+plot/populationVectorCorrelation.m b/+plot/populationVectorCorrelation.m new file mode 100644 index 0000000..340dae7 --- /dev/null +++ b/+plot/populationVectorCorrelation.m @@ -0,0 +1,60 @@ + + + +for k =1:length(areas) + meanSignal_area(k).deconv = []; + for j = 1:length(data(areas(k)).sessionIDs) + for f = 1:length(mData(areas(k),j).LNP_results) + if ~isempty(find(mData(areas(k),j).LNP_results(f).selected_model == 1)) && isempty(find(mData(areas(k),j).LNP_results(f).selected_model == 5)) + meanSignal_area(k).deconv = [ meanSignal_area(k).deconv; rmaps(areas(k),j).deconv(f,:)]; + end + end + end + + + for i = 1:size(meanSignal_area(k).deconv,2) + for j =1:size(meanSignal_area(k).deconv,2) + f = corrcoef(meanSignal_area(k).deconv(:,i), meanSignal_area(k).deconv(:,j),'rows','complete'); + corrval(i,j) = f(1,2); + end + end + + + fig = figure(); + imagesc(flipud(corrval)) + colormap('jet') + limsY = get(gca,'YLim'); + hold on + for i = 1: length(cue1) + plot(gca, [cue1(i) cue1(i)],limsY, 'LineStyle','--','Color', [1 1 1], 'LineWidth', 2); + end + + for i = 1: length(cue2) + plot(gca, [cue2(i) cue2(i)],limsY, 'LineStyle','--','Color', [1 1 1], 'LineWidth', 2); + end + limsX = get(gca,'XLim'); + + for i = 1: length(cue1) + plot(gca,limsX, [limsX(2)-cue1(i) limsX(2)-cue1(i)], 'LineStyle','--','Color', [1 1 1], 'LineWidth', 2); + end + + for i = 1: length(cue2) + plot(gca,limsX, [limsX(2)-cue2(i) limsX(2)-cue2(i)], 'LineStyle','--','Color', [1 1 1], 'LineWidth', 2); + end + % + xlabel('position (cm)'); ylabel('position (cm)'); + % % + xticks([1 20 40 60 78]); + xticklabels({'0','40','80','120', '157'}); + set(gca, 'FontSize',18, 'FontName', 'Arial') + % ylabel('# axons') + % + yticks(flipud([78-1 78-20 78-40 78-60 1]')); + yticklabels({'157','120','80','40', '0'}'); + title(data(areas(k)).area) + colorbar + +end + + + diff --git a/+plot/traces.m b/+plot/traces.m new file mode 100644 index 0000000..1158887 --- /dev/null +++ b/+plot/traces.m @@ -0,0 +1,55 @@ + + +sData = mData(area,session).sData; +figTraces = figure(); +set(figTraces, 'Units', 'centimeters', 'Position', [0 0 20 30]); + +ax_signal = axes('Position',[0.1 0.31 0.8 0.65]); hold(ax_signal, 'on') +ax_pos = axes('Position',[0.1 0.19 0.8 0.1]); +ax_speed = axes('Position',[0.1 0.07 0.8 0.1]); + + +% fig = gcf; + +nRois = length(rois); +signal = sData.imdata.roiSignals(2).dff; +signalDeconv = sData.imdata.roiSignals(2).deconv; + +added = 1; + +% adjustSizeDFF = 1200; +% adjustSizeDeconv = 3500; + +spaceBetweenCells = 1; + +title(data(area).area) + +jj = 1; minVal = []; maxVal = []; +for rNo = 1:nRois + signalTemp = signal(rois(rNo),:)-spaceBetweenCells*(rNo); + signalDeconvTemp = signalDeconv(rois(rNo),:)-spaceBetweenCells*(rNo)-0.1*spaceBetweenCells; + + plot(ax_signal, signalTemp, 'LineWidth',0.5,'Color',[0 0 0]); jj = jj+1; + plot(ax_signal, signalDeconvTemp,'Color','r'); + + maxVal = max([maxVal max(signalTemp)]); + minVal = min([minVal min(signalDeconvTemp)]); + jj = jj+1; +end + +ax_signal.XAxis.Visible = 'off'; ax_signal.YAxis.Visible = 'off'; +plot(ax_pos,sData.behavior.wheelPosDs,'k', 'LineWidth',1); +plot(ax_speed,sData.behavior.runSpeedDs,'k', 'LineWidth',1); + + +set(ax_pos, 'FontSize',12, 'FontName', 'Arial','YTickLabels',[],'XTickLabels',[], 'Box', 'off', 'XLim', xlimit); +set(ax_speed, 'FontSize',12, 'FontName', 'Arial','YTickLabels',[],'XTickLabels',[], 'Box', 'off', 'XLim', xlimit); +set(ax_signal, 'FontSize',12, 'FontName', 'Arial','YTickLabels',[],'XTickLabels',[], 'Box', 'off', 'XLim', xlimit); + +ax_pos.YLabel.String= 'pos'; +ax_speed.YLabel.String= 'speed'; + + +ax_speed.XTick = xlimit(1):31*5:xlimit(2); +ax_speed.XTickLabel = 0:5:5*length(ax_speed.XTick); +ax_speed.XLabel.String = 'Time(s)'; \ No newline at end of file diff --git a/Figures/create_all_figures.m b/Figures/create_all_figures.m new file mode 100644 index 0000000..f73a49d --- /dev/null +++ b/Figures/create_all_figures.m @@ -0,0 +1,11 @@ + +load_all_data() % find explanation of sData structure in: sData_template.m +load_LNP_results() % if not available run analysis/cellClassification +load_decoding_results() % if not available run decoding +load_skaggs_results() % if not available run skaggs +helper.createRMaps; + +create_figure_1; +create_figure_2; +create_figure_3; + diff --git a/Figures/create_figure_1.m b/Figures/create_figure_1.m new file mode 100644 index 0000000..8f5979b --- /dev/null +++ b/Figures/create_figure_1.m @@ -0,0 +1,8 @@ + + +fig_1_b_c; +fig_1_f; +fig_1_g_o; +fig_1_p; +fig_1_l; +fig_1_h_n; \ No newline at end of file diff --git a/Figures/create_figure_2.m b/Figures/create_figure_2.m new file mode 100644 index 0000000..aafe395 --- /dev/null +++ b/Figures/create_figure_2.m @@ -0,0 +1,11 @@ + + +fig_2_b; +fig_2_c; +fig_2_d_e; +fig_2_g; +fig_2_h; + + + + diff --git a/Figures/create_figure_3.m b/Figures/create_figure_3.m new file mode 100644 index 0000000..de0f2f6 --- /dev/null +++ b/Figures/create_figure_3.m @@ -0,0 +1,6 @@ + +fig_3_a; +fig_3_b_c_i; +fig_3_e_f; +fig_3_d_g; +fig_3_h; \ No newline at end of file diff --git a/Figures/figure_1/fig_1_b_c.m b/Figures/figure_1/fig_1_b_c.m new file mode 100644 index 0000000..ab7df42 --- /dev/null +++ b/Figures/figure_1/fig_1_b_c.m @@ -0,0 +1,80 @@ + +k = 1; +for i = 1:length(data) + for f = 1:length(data(i).sessionIDs ) + + sData = mData(i,f).sData; + + % find mean lick frequency across position + [lickFrequencyDs,lickFrequency] = helper.extractLickFrequency(sData, 10000, 15, 15); + + matrixLick = plot.createHeatMap(sData.trials.trialLength, sData.behavior.wheelPosDsBinned,lickFrequencyDs); + meanLick{i}(f,:) = nanmean(matrixLick); + + % find first licks + firstLicks{k} = []; + for r = 1:size(matrixLick,1) + idxLick = find(matrixLick(r,15:end)>1); + if ~isempty(idxLick) + firstLicks{k} = [firstLicks{k}, 15+idxLick(1)]; + end + end + k = k+1; + + + %% create run matrix + + matrixRun = plot.createHeatMap(sData.trials.trialLength, sData.behavior.wheelPosDsBinned,sData.behavior.runSpeedDs); + meanRun{i}(f,:) = nanmean(matrixRun); + + end +end + + + + + F = 1:size(meanRun{1},2); + +c1 = [80 134 198]/255; +c2 = [0 0 0]/255; + +c = [linspace(c1(1),c2(1),2)', linspace(c1(2),c2(2),2)', linspace(c1(3),c2(3),2)']; + +%% plot run behavior +fig_1_a = figure(); +amean = nanmean(cell2mat(meanRun')); +astd = nanstd(cell2mat(meanRun')); +plot(amean,'LineWidth',1,'Color',c(1,:)) +hold on +patch([F(~isnan(amean)) fliplr(F(~isnan(amean)))],[amean(~isnan(amean))+astd(~isnan(amean)) fliplr(amean(~isnan(amean))-astd(~isnan(amean)))],c(1,:),'linestyle','none','FaceAlpha', 0.2); + +box off +set(gca,'FontName', 'Arial', 'FontSize', 12) +xlabel(gca,'Position (cm)') +ylabel(gca,'Speed (cm/s)') +xticks(gca,[1 25 50 78]) +xticklabels(gca,{'0' ,'50','100','157'}) +xlim([1 79]) + +%% plot lick behavior +fig_1_b = figure(); +amean = nanmean(cell2mat(meanLick')); +astd = nanstd(cell2mat(meanLick'))./sqrt(8); +plot(amean,'LineWidth',1,'Color',c(2,:)) +hold on +patch([F(~isnan(amean)) fliplr(F(~isnan(amean)))],[amean(~isnan(amean))+astd(~isnan(amean)) fliplr(amean(~isnan(amean))-astd(~isnan(amean)))],c(2,:),'linestyle','none','FaceAlpha', 0.2); + +firstLicksVal = []; +for i = 1:length(firstLicks); firstLicksVal = [firstLicksVal, nanmean(firstLicks{i})]; end + +plot(gca, nanmean(firstLicksVal),6,'.','MarkerSize', 25,'Color', 'k') +errorbar(gca,nanmean(firstLicksVal),6,nanstd(firstLicksVal),'horizontal','k','LineWidth',1.5, 'Color', c(2,:)) + +box(gca, 'off') +set(gca, 'FontName', 'Arial', 'FontSize', 12) +xlabel(gca,'Position (cm)') +ylabel(gca,'Lick rate (Hz)') +xticks(gca,[1 25 50 78]) +xticklabels(gca,{'0' ,'50','100','157'}) +xlim([1 79]) + diff --git a/Figures/figure_1/fig_1_f.m b/Figures/figure_1/fig_1_f.m new file mode 100644 index 0000000..875e801 --- /dev/null +++ b/Figures/figure_1/fig_1_f.m @@ -0,0 +1,8 @@ + + +area = 7; +session = 9; +rois = [40 82 92]; +xlimit = [1000 2000]; + +plot.traces; \ No newline at end of file diff --git a/Figures/figure_1/fig_1_g_o.m b/Figures/figure_1/fig_1_g_o.m new file mode 100644 index 0000000..364617e --- /dev/null +++ b/Figures/figure_1/fig_1_g_o.m @@ -0,0 +1,5 @@ + +% requires mData, data and rmaps +% helper.createRMaps; +areas = [7 2]; +plot.plotSortedSequence; \ No newline at end of file diff --git a/Figures/figure_1/fig_1_h_n.m b/Figures/figure_1/fig_1_h_n.m new file mode 100644 index 0000000..958887b --- /dev/null +++ b/Figures/figure_1/fig_1_h_n.m @@ -0,0 +1,5 @@ + +areas = [7,2]; + +plot.fraction_tuning; + diff --git a/Figures/figure_1/fig_1_i.m b/Figures/figure_1/fig_1_i.m new file mode 100644 index 0000000..f849503 --- /dev/null +++ b/Figures/figure_1/fig_1_i.m @@ -0,0 +1,19 @@ + +% initiate variables + +figure() +area = 7; +plot.decodingVsNoCells + +% plot results + +plot(gca,cellNoSum,nanmean(rmse_all,1), 'Color',[0 0 0], 'LineWidth', 1.5) +yval = rmse_all; +patch(gca,[cellNoSum fliplr(cellNoSum)],[nanmean(yval)+nanstd(yval)./sqrt(size(rmse_all,1)) fliplr(nanmean(yval)-nanstd(yval)./sqrt(size(rmse_all,1)))],[0 0 0],'linestyle','none','FaceAlpha', 0.2); + +ylabel('Decoding error (cm)') +xlabel('# Cells') +set(gca,'FontName','Arial','FontSize',12) +box off +ylim([0 55]) + diff --git a/Figures/figure_1/fig_1_l.m b/Figures/figure_1/fig_1_l.m new file mode 100644 index 0000000..96df443 --- /dev/null +++ b/Figures/figure_1/fig_1_l.m @@ -0,0 +1,11 @@ +% M2 +% sData = load('/Users/annachristinagarvert/UIO Physiology Dropbox Dropbox/Lab Data/Michele Gianatti/Data&Analysis/sData_SNR/M2/m5094-20200522-03/final_sData.mat'); +% rois = [234 243 254]; +% xlimit = [4500 5400]; + +area = 2; +session = 3; +rois = [234 243 254]; +xlimit = [800 2000]; + +plot.traces; \ No newline at end of file diff --git a/Figures/figure_1/fig_1_m.m b/Figures/figure_1/fig_1_m.m new file mode 100644 index 0000000..37fef15 --- /dev/null +++ b/Figures/figure_1/fig_1_m.m @@ -0,0 +1,51 @@ + +area = 2; % M2 +session = 3; +roi = 243; +% isPlaceCell = mData(area,session).LNP_results(roi); check model performance + +signal =normalize(mData(area,session).sData.imdata.roiSignals(2).dff(roi,:),'range',[0,1]); +binnedMatrix = plot.createHeatMap(mData(area,session).sData.trials.trialLength,mData(area,session).sData.behavior.wheelPosDsBinned,signal); + +fig_1_m = figure(); +imagesc(binnedMatrix) +set(gca, 'FontName', 'Arial', 'FontSize', 12) +xticks(gca,[1 20 40 60 78]); +xticklabels(gca,{'0', '40','80','120','158'}); +colormap(jet) +ylabel('Trials'); +xlabel('Position (cm)'); +colormap(jet); +xlim([1 78]) +title('Normalized DFF') +colorbar() + + + + +function matrix = createHeatMap(lapVector,position,behavior) + +heatmap = zeros(max(lapVector)+1,ceil(max(position))+1); +occupancymap = zeros(max(lapVector)+1,ceil(max(position))+1); +offset = min(lapVector); +for x = 1:length(position) + x_pos = round(position(x)); + y_pos = lapVector(x)+1+(-offset); + heatmap(y_pos,x_pos+1) = heatmap(y_pos,x_pos+1) + behavior(x); + occupancymap(y_pos,x_pos+1) = occupancymap(y_pos,x_pos+1) + 1; +end + +matrix = heatmap./occupancymap; + + +% sometimes, the positions are not evenly sampled and some values in +% the matrix remain not assigned (here i fill in these with +% neighbourign values) +% matrix = matrix'; +% while ~isempty(find(isnan(matrix))) +% matrix(isnan(matrix)) = matrix(find(isnan(matrix))-1); +% end +% matrix = matrix'; + + +end \ No newline at end of file diff --git a/Figures/figure_1/fig_1_p.m b/Figures/figure_1/fig_1_p.m new file mode 100644 index 0000000..4a9d0f6 --- /dev/null +++ b/Figures/figure_1/fig_1_p.m @@ -0,0 +1,3 @@ +areas = 2; + +plot.populationVectorCorrelation \ No newline at end of file diff --git a/Figures/figure_2/fig_2_b.m b/Figures/figure_2/fig_2_b.m new file mode 100644 index 0000000..ae8e9a0 --- /dev/null +++ b/Figures/figure_2/fig_2_b.m @@ -0,0 +1,46 @@ + +% OFC +area = 1; +session = 5; +rois = [12 13 15]; +xlimit = [12160 13160]; + +plot.traces; + + +% ACC +area = 3; +session = 4; +rois = [145 165 175]; +xlimit = [7000 8500]; + +plot.traces; + + +% Thalamus +area = 4; +session = 16; +rois = [43 48 270]; +xlimit = [4000 6000]; + +plot.traces; + +% PPC +area = 5; +session = 1; +rois = [204 216 217]; +xlimit = [5200 8000]; + +plot.traces; + + +% V2 + +area = 6; +session = 1; +rois = [274 164 329]; +xlimit = [2500 5000]; + +plot.traces; + + diff --git a/Figures/figure_2/fig_2_c.m b/Figures/figure_2/fig_2_c.m new file mode 100644 index 0000000..5203d23 --- /dev/null +++ b/Figures/figure_2/fig_2_c.m @@ -0,0 +1,8 @@ + +% requires mData, data and rmaps +% helper.createRMaps; make sure this ran before, this crease rmaps +areas = [1 3 4 5 6]; +plot.plotSortedSequence; + + + diff --git a/Figures/figure_2/fig_2_d_e.m b/Figures/figure_2/fig_2_d_e.m new file mode 100644 index 0000000..ed7a712 --- /dev/null +++ b/Figures/figure_2/fig_2_d_e.m @@ -0,0 +1,70 @@ + + + + +%% plot fraction of place tuned axons +fig_2_d = figure(); + +sorting = [3 5 4 2 6 1]; % this sorts the areas with increasing place tuned fraction +for l = 1:6 + for i = 1: length(data(l).sessionIDs) + pcNo = 0; + for k = 1:length(mData(l,i).LNP_results) + if ~isempty(find(mData(l,i).LNP_results(k).selected_model == 1, 1)) + pcNo = pcNo +1; + end + end + + pc{l}(i) = 100*pcNo/length(mData(l,i).LNP_results); + scatter(sorting(l), pc{l}(i), 70,[0.5 0.5 0.5]);%col(l,:)) + hold on + end + PCall(l) = nanmean(pc{l}); + PCallErr(l) = nanstd(pc{l})/sqrt(length(pc{l})); + + scatter(sorting(l),PCall(l),70,[0 0 0],'filled') + errorbar(sorting(l),PCall(l),PCallErr(l),'k','LineWidth',1.5) +end + +xticks([1:6]) +xlim([0.5 6.5]) +set(gca,'FontName','Arial','FontSize',12) +xticklabels({'V2','Thal','OFC','ACC','M2','PPC'}) +ylabel('Position tuned axon (%)') + +%% check p values figure 2 e +sorting_inverse = [6 4 1 3 2 5]; + +for i = 1:6 + for j = 1:6 + pval(i,j) = ranksum(pc{sorting_inverse(i)},pc{sorting_inverse(j)}); + end +end + +matr =triu(ones(6,6)); +matr(eye(6,6)==1) = NaN; +matr(matr== 0)= NaN; +pval(isnan(matr)) = NaN; +[corrected_p, h]=helper.bonf_holm(pval); + + +fig_2_e = figure(); +imagesc(triu(corrected_p)) +hold on +alpha(matr) + +xticklabels({'V2','Thal','OFC','ACC','M2', 'PPC'}) +yticklabels({'V2','Thal','OFC','ACC','M2', 'PPC'}) + + +c = colorbar('Limits',[0 0.05]); +caxis([0 0.05]) +ax = gca; +ax.FontSize = 12; +ax.FontName = 'Arial'; +ax.XColor = [0 0 0]; +ax.YColor = [0 0 0]; +ax.XLabel.Color = [0 0 0]; +ax.YLabel.Color = [0 0 0]; + + \ No newline at end of file diff --git a/Figures/figure_2/fig_2_g.m b/Figures/figure_2/fig_2_g.m new file mode 100644 index 0000000..5bf7c25 --- /dev/null +++ b/Figures/figure_2/fig_2_g.m @@ -0,0 +1,62 @@ + +% helper.createRMaps; + +meanSignal_area = struct(); + +for k =1:length(data) + + meanSignal_area(k).dff = []; + + for j = 1:length(data(k).sessionIDs) + for f = 1:length(mData(k,j).LNP_results) + if ~isempty(find(mData(k,j).LNP_results(f).selected_model == 1)) && isempty(find(mData(k,j).LNP_results(f).selected_model == 5)) + meanSignal_area(k).dff = [ meanSignal_area(k).dff; rmaps(k,j).dff(f,:)]; + end + end + + end +end + + +col = [122 172 210 + 80 106 139 + 217 85 88 + 129 89 162 + 180 151 94 + 179 205 142 + 0 0 0 + 0 0 0]./255; + + +fig_g = figure(); +for k = 1:length(area)-2 + + [~, idx] = max(meanSignal_area(k).dff, [], 2); % Returns 2 arrays. Val is the max value of each ROI, idx is the column index of the max value of each ROI. + idxMatrix = horzcat(idx, meanSignal_area(k).dff); % Appends idx column array to the ROI matrix. + sortedMatrix = sortrows(idxMatrix, 1); % Sorts ROI matrix by the index of their max values in ascending order. + sortedMatrix(: , 1) = []; % Removes indexes from first column of sortedMatrix. + + sortedMatrix_aligned = NaN(size(sortedMatrix)); + + for m = 1:length(idx) + sortedMatrix_aligned(m,:) = circshift(meanSignal_area(k).dff(m,:), 78-idx(m)+39); + end + + plot(nanmean(sortedMatrix_aligned), 'LineWidth', 1.5, 'Color', col(k,:));hold on + +end +ax.FontName = 'Arial'; +ax.XColor = [0 0 0]; +ax.YColor = [0 0 0]; +ax.XLabel.Color = [0 0 0]; +ax.YLabel.Color = [0 0 0]; +box off +legend({'OFC', 'M2', 'ACC', 'Thal', 'PPC', 'V2'}) +xlabel('Position relative to reward(cm)'); +% +xticks([1 19 39 58 78]); +xticklabels({'-80','-40','0','40','78'}); +set(gca,'FontSize',12, 'FontName', 'Arial') +ylabel('Deconv DFF') +hold on + diff --git a/Figures/figure_2/fig_2_h.m b/Figures/figure_2/fig_2_h.m new file mode 100644 index 0000000..249ed24 --- /dev/null +++ b/Figures/figure_2/fig_2_h.m @@ -0,0 +1,83 @@ + + +for l = 1:length(data)-2 + placeCellSignal = []; %width{l} = []; no_PF{l} = []; + for i = 1: length(data(l).sessionIDs) + width{l}{i} = []; + no_PF{l}{i} = []; + sData = mData(l,i).sData; + + pcNo = []; + for k = 1:length(mData(l,i).LNP_results) + if ~isempty(find(mData(l,i).LNP_results(k).selected_model == 1)) && isempty(find(mData(l,i).LNP_results(k).selected_model == 5)) + pcNo = [pcNo k]; + end + end + + if ~isempty(pcNo) + for r = 1:length(pcNo) + signal_matrix = helper.splitInTrialsPos(sData.imdata.roiSignals(2).dff, sData,pcNo(r)); + tuning_curve = nanmean(signal_matrix,1); + + halfMaximum = min(tuning_curve)+(max(tuning_curve)-min(tuning_curve))/2; + + idxVals = find(tuning_curve > halfMaximum); + widthTemp = []; + + m = 1; f = 1; widthTemp{f} = 1; tempM = idxVals(1); + while m < length(idxVals) + if (idxVals(m+1) - idxVals(m))>1 + widthTemp{f} = idxVals(m)-tempM+1; + tempM = idxVals(m+1); + f = f+1; + end + m = m+1; + end + widthTemp{f} = idxVals(m)-tempM+1; + + % correct for circularity + if length(widthTemp) >= 2 + if ~isempty(find(idxVals == 1)) && ~isempty(find(idxVals == 78)) + widthTemp{1} = widthTemp{1}+widthTemp{end}; + widthTemp(end) = []; + end + end + + widthTemp2 =cell2mat(widthTemp); + + widthTemp2(widthTemp2 < 5) = []; + width{l}{i} = [width{l}{i} widthTemp2]; + no_PF{l}{i}(end+1) = length(widthTemp2); + + end + end + + end +end + + + +for i = 1:length(data)-2 + for k = 1:length(width{i}) + widthMean{i}(k) = nanmean(width{i}{k}); + end + widthMeanArea(i) = nanmean(widthMean{i}); + widthSEDArea(i) = nanstd(2*widthMean{i})/sqrt(length(width{i})); +end + + +fig_h = figure(); +for i = 1:length(data)-2 + scatter(i,2* widthMeanArea(i),80, 'k', 'filled') ;hold on + errorbar(i, 2*widthMeanArea(i),widthSEDArea(i),'LineWidth',1.5, 'Color','k') + +end + +xticklabels(gca,{'OFC', 'M2', 'ACC', 'Thal','PPC', 'V2'}) +xticks(gca,[1:6]) +xlim(gca,[0 7]) +ylabel('Place field width(cm)') +ylim(gca,[0 100]) +set(gca,'FontName','Arial','FontSize',12) + + \ No newline at end of file diff --git a/Figures/figure_3/fig3_b_c_i.m b/Figures/figure_3/fig3_b_c_i.m new file mode 100644 index 0000000..ebd7d9e --- /dev/null +++ b/Figures/figure_3/fig3_b_c_i.m @@ -0,0 +1,159 @@ + + +col = [122 172 210 + 80 106 139 + 217 85 88 + 129 89 162 + 180 151 94 + 179 205 142 + 0 0 0 + 255 50 80]./255; + +randAll = []; + +type ='d_data'; + +maxval =154.3735; +classes = 2.6165:5.233:157.05; +classbin = 1:30; + +rmse_vs_pos = cell(6,1); + +fig_i = figure(); +for r = 1:6 + for f = 1:length(data(r).sessionIDs) + rmse_vs_pos_folds = NaN(5,length(classes)); + pred_PosMean = NaN(5,length(classes)); + + if ~isempty(mData(r,f).(type)) + for l = 1:5 + pred_pos = mData(r,f).(type)(length(mData(r,f).(type))).iter{1,l}.predPos_test; + actual_class = mData(r,f).(type)(length(mData(r,f).(type))).iter{1, l}.realPos_test; + + % transform back from class to position bin + actual_pos = NaN(length(actual_class),1); + for t = 1:length(actual_class) + [~,ind] = min(abs(classes-actual_class(t))); + actual_pos(t) = classes(ind); + end + + % calculate error for all time bins + predError = zeros(length(pred_pos),1); + for i = 1:length(actual_pos) + predError(i) = abs(actual_pos(i)-pred_pos(i)); + if predError(i) > maxval/2 && actual_pos(i) > pred_pos(i) + predError(i) = abs(maxval-actual_pos(i)+pred_pos(i)); + pred_pos(i) = maxval+pred_pos(i); + elseif predError(i) > maxval/2 &&actual_pos(i) < pred_pos(i) + predError(i) = abs(maxval-pred_pos(i)+actual_pos(i)); + pred_pos(i) = -(maxval-pred_pos(i)); + end + end + + + % calculate error for each class (position related error + for k= 1:length(classes) + idx = find(actual_pos == classes(k)); + rmse_vs_pos_folds(l,k) = nanmean(predError(idx')); + end + + + end + + + rmse_vs_pos{r}(f,:) = nanmean(rmse_vs_pos_folds); + rmse(r,f) = nanmean(predError); + + end + end + + + % plot rmse vs pos + F = 1:31; + amean = nanmean(rmse_vs_pos{r},1); + amean = [amean amean(1)]; + astd = nanstd(rmse_vs_pos{r})./sqrt(size(rmse_vs_pos{r},1)); + astd = [astd astd(1)]; + patch(gca,[F(~isnan(amean)) fliplr(F(~isnan(amean)))],[amean(~isnan(amean))+astd(~isnan(amean)) fliplr(amean(~isnan(amean))-astd(~isnan(amean)))],col(r,:),'linestyle','none','FaceAlpha', 0.2); + hold on + plot(amean, 'LineWidth', 1.2, 'Color', col(r,:)) + +end + + +ylabel('Decoding error (cm)') +xlabel('Position (cm)') +set(gca, 'FontName', 'Arial', 'FontSize', 12) +box(gca, 'off') +xticks(gca, [1 7.5 15 22.5 30]) +xticklabels(gca, {'0', '40','80' ,'120','157'}) + + +cues = [42 47 122 128 63 68 100 106]./5.2333; + +for i = 1:8 + xline(cues(i),'--',[0.5 0.5 0.5],'LineWidth',1.5) +end + + +% plot mean rmse for each area +fig_b_or_e = figure(); +sorting = [6 1 4 2 3 5]; +inverse_sorting = [2 4 5 3 6 1]; + +sorted_x_labels = cell(6,1); + +for r = 1:6 + amean = nanmean(rmse_vs_pos{r},2); + scatter(ones(1,length(amean))*sorting(r), amean, 70,[0.5 0.5 0.5]);%col(l,:)) + hold( 'on') + + scatter(sorting(r), nanmean(amean), 70,[0 0 0],'filled');%col(l,:)) + errorbar(sorting(r), nanmean(amean),nanstd(amean)/sqrt(length(amean)),'Color' ,'k', 'LineWidth', 1.5) + sorted_x_labels{r} = data(inverse_sorting(r)).area; +end + +xlim([0.5 6.5]) +xticks(gca, 1:6) +xticklabels(gca, sorted_x_labels) +xtickangle(gca, 45) +set(gca, 'FontName', 'Arial', 'FontSize', 12) +ylabel(gca,'Decoding error (cm)') +box(gca, 'off') + + + +% check wether rmse are singificantly different between areas +% - bonferroni holm corrected pvalues +rmse(rmse==0) = NaN; + +for i = 1:6 + for j = 1:6 + pval(i,j) = ranksum(rmse(inverse_sorting(i),:), rmse(inverse_sorting(j),:)); + end + +end + +matr =triu(ones(6,6)); +matr(eye(6,6)==1) = NaN; +matr(matr== 0)= NaN; +pval(isnan(matr)) = NaN; +[corrected_p, h]=helper.bonf_holm(pval); +imagesc(triu(corrected_p)) +hold on +alpha(matr) + +c = colorbar('Limits',[0 0.05]); +caxis([0 0.05]) +ax = gca; +ax.FontSize = 12; +ax.FontName = 'Arial'; +ax.XColor = [0 0 0]; +ax.YColor = [0 0 0]; +ax.XLabel.Color = [0 0 0]; +ax.YLabel.Color = [0 0 0]; + +xticklabels(sorted_x_labels); +yticklabels(sorted_x_labels); + + diff --git a/Figures/figure_3/fig3_d_g_old.m b/Figures/figure_3/fig3_d_g_old.m new file mode 100644 index 0000000..27e9af7 --- /dev/null +++ b/Figures/figure_3/fig3_d_g_old.m @@ -0,0 +1,579 @@ +% close all +% clear all + + +% MP.figure(4) +% axes = findall(gcf,'Type','Axes'); +%direction = '/Users/annachristinagarvert/Desktop/Data/Analysed 03:20/sData/'; +%% LDDM m5051 #4 m5088 #7 m5074 #4 m5052 #1 m5075 #4 m5087 #5 + +% files{7} = {'m5115-20201203-02','m5115-20201203-01','m5115-20201202-01','m5114-20201203-01',... +% 'm5114-20201202-02','m5114-20201202-01','m5108-20201120-03','m5108-20201120-02','m5108-20201120-01'}; + +% files{1} = {'m5142-20210709-01','m5143-20210710-02'}; +% +% %% M2 m5094, m5095 +% % could be added: 'm5095-20200520-01_M2' +% files{2} = {'m5094-20200520-01','m5095-20200520-01','m5094-20200522-03',... +% 'm5095-20200522-01','m5128-20210324-01','m5122-20210304-01',... +% 'm5122-20210305-01'}; +% +% +% files{3} = {'m5082-20200312-01','m5082-20200312-02',... +% 'm5085-20200313-02','m5085-20200314-01','m5124-20210323-02','m5124-20210324-02',... +% 'm5125-20210324-01','m5125-20210323-01', 'm5125-20210325-01'}; +% %m5082-20200313-02 'm5085-20200313-01' +% +% files{4} = {'m5074-20191202-01','m5075-20191203-01', 'm5074-20191205-01','m5074-20191219-01',... +% 'm5075-20191205-01','m5075-20191205-02','m5075-20191219-01','m5088-20200430-02','m5088-20200430-03',... +% 'm5088-20200504-01','m5088-20200504-02','m5088-20200504-03','m5088-20200508-01','m5088-20200508-03',... +% 'm5140-20210709-01','m5140-20210710-01','m5140-20210710-02','m5141-20210709-02','m5141-20210710-01'}; +% % +% files{5} = {'m5121-20210204-03','m5121-20210205-01','m5121-20210205-02','m5119-20210204-01',... +% 'm5119-20210204-02','m5153-20211208-01',... +% 'm5153-20211208-02','m5153-20211209-02'}; +% +% +% files{6} = {'m5137-20210610-01','m5137-20210610-02',... +% 'm5138-20210610-01', 'm5138-20210610-02','m5138-20210611-01',... +% 'm5139-20210611-01','m5139-20210611-02','m5139-20210612-01',... +% 'm5146-20211021-01','m5146-20211024-02','m5146-20211022-01',... +% 'm5147-20211022-01','m5147-20211024-01'}; +% +% files{7} = {'m5108-20201120-01','m5108-20201120-02','m5108-20201120-03','m5114-20201202-01','m5114-20201202-02',... +% 'm5114-20201203-01','m5115-20201202-01','m5115-20201203-01','m5115-20201203-02'}; +% +% % files{8} = {'m5091-20200504-03','m5091-20200508-03','m5091-20200508-06',... +% % 'm5102-20200713-01'}; + +% load_all_data() + + + +col = [122 172 210 + 80 106 139 + 217 85 88 + 129 89 162 + 180 151 94 + 179 205 142 + 0 0 0 + 0 0 0]./255; % OFC, M2, ACC, Thal, PPC, V2 +% % LDDM, ACC, M2, PPC, V2, OFC +% Tha, ACC, M2, PPC, V2, OFC +% col = [0 0 0; 129 89 162 ; 217 85 88; 80 106 139; 180 151 94; 179 205 142; 122 172 210]./255; + +% M2 [238 155 56]./255 +% RSC [235 72 83]./255 +% LDDM [77 52 142]./255 +% ACC [119 152 107]./255 +%% separate plots for mice: +% figMeanAll = figure; +% fig = figure(); +% hold(axes(4), 'off') +% area = {'OFC','M2','ACC','Thal', 'PPC', 'V2','RSC','V1'}; +area = {'OFC','M2','ACC','Thal', 'PPC', 'V2','RSC'}; + +% figure() + +% remove_OFC_sessions= [6,10,11,12]; +% file{1}(remove_OFC_sessions)= []; + +for r = 1:7 %length(area) + rmse_area{r} = []; lengthCells = [];cellsSpaced = []; rand_pos{r}= []; rmse_rand{r}= [];rmse_rand_CI{r}= []; + cellNo{r} = []; + f = 1; + for l = 1:length(file{r}) + + load(fullfile('/Users/annachristinagarvert/Desktop/Michele/analysis/LNP/5Folds/',area{r},file{r}{l},'axon_final.mat')) + pcNo = 0; + for p = 1:length(axon) + if ~isempty(find(axon(p).selected_model == 1)) + pcNo = pcNo +1; + end + end + + % there has toobe a least one pc + if pcNo > 0 + +% load(strcat('/Users/annachristinagarvert/Desktop/Michele/analysis/decoding/decodingVsNoCells/',area{r},'/',file,'/d_data_varN.mat')) + + load(strcat('/Users/annachristinagarvert/UIO Physiology Dropbox Dropbox/Lab Data/Michele Gianatti/Data&Analysis/Paper/decoding/decodingOnlyPC/',area{r},'/',file{r}{l},'/d_data_varN.mat')) + load(fullfile('/Users/annachristinagarvert/UIO Physiology Dropbox Dropbox/Lab Data/Michele Gianatti/Data&Analysis/Paper/sData/',area{r},'/',file{r}{l},'/final_sData.mat')) + + + predPos = d_data_varN(end).iter{1, d_data_varN(end).medianTestFold}.predPos_test; + realPos = d_data_varN(end).iter{1, d_data_varN(end).medianTestFold}.realPos_test; + + maxval =154.3735;% max([max(predPos(:)) max(realPos(:))]); + predError = zeros(length(realPos),1); + for k = 1:length(realPos) + predError(k) = abs(realPos(k)-predPos(k)); + if predError(k) > maxval/2 && realPos(k) > predPos(k) + predError(k) = abs(maxval-realPos(k)+predPos(k)); + elseif predError(k) > maxval/2 && realPos(k) < predPos(k) + predError(k) = abs(maxval-predPos(k)+realPos(k)); + end + end + + meanRMSE(r,f) = nanmean(predError); + + % for m = 1:5000 + % + % randPos = circshift(predPos, randperm(length(predPos),1)); + % + % maxval =154.3735;% + % randError = zeros(length(realPos),1); + % for k = 1:length(realPos) + % randError(k) = abs(realPos(k)-randPos(k)); + % if randError(k) > maxval/2 && realPos(k) > randPos(k) + % randError(k) = abs(maxval-realPos(k)+randPos(k)); + % elseif randError(k) > maxval/2 && realPos(k) < randPos(k) + % randError(k) = abs(maxval-randPos(k)+realPos(k)); + % end + % end + % + % rand_rmse(m) = nanmean(randError); + % + % end + % + % mean_rand(r,f) = nanmean(rand_rmse); + % tempCI = quantile(rand_rmse,[0.05 0.95]); + % mean_rand_CI(r,f) = tempCI(1); + f = f+1; + end + end + +end + +meanRMSE(meanRMSE==0) = NaN; +% mean_rand(mean_rand == 0) = NaN; +% mean_rand_CI(mean_rand_CI == 0) = NaN; + +figure() +dir = [7 2 4 5 3 6 1]; % only PC xticklabels({'RSC','M2','PPC', 'ACC','Thal', 'V2','OFC'}) + +% dir = [7 2 5 4 3 6 1]; +for r = 1:6 + for l = 1:size(meanRMSE,2) + scatter(dir(r),meanRMSE(r,l),70, [0.5 0.5 0.5]); + hold on + end + scatter(dir(r),nanmean(meanRMSE(r,:)),70, [0 0 0],'filled'); + errorbar(dir(r),nanmean(meanRMSE(r,:)),nanstd(meanRMSE(r,:))/sqrt(sum(~isnan(meanRMSE(r,:)))),'k','LineWidth',1.5); +end +% xticklabels({'RSC','M2','PPC', 'Thal','ACC', 'V2','OFC'}) + +xlim([1 8]) +ylim([5 55]) +xticks([2:7]) +ax = gca; +ax.FontSize = 12; +ax.FontName = 'Arial'; +ax.XColor = [0 0 0]; +ax.YColor = [0 0 0]; +ax.XLabel.Color = [0 0 0]; +ax.YLabel.Color = [0 0 0]; +xticklabels({'M2','PPC','ACC', 'Thal', 'V2','OFC'}) +% yline(nanmean(mean_rand_CI(:))) +% yline(nanmean(mean_rand(:))) + +sortedRMSE{1} = meanRMSE(2,:); +sortedRMSE{2} = meanRMSE(5,:); +sortedRMSE{3} = meanRMSE(3,:); +sortedRMSE{4} = meanRMSE(4,:); +sortedRMSE{5} = meanRMSE(6,:); +sortedRMSE{6} = meanRMSE(1,:); + +for i = 1:6 + for j = 1:6 + pval(i,j) = ranksum(sortedRMSE{i},sortedRMSE{j}); + end +end + +matr =triu(ones(6,6)); +matr(eye(6,6)==1) = NaN; +matr(matr== 0)= NaN; +pval(isnan(matr)) = NaN; +[corrected_p, h]=bonf_holm(pval); + +figure() +imagesc(triu(corrected_p)) +hold on +alpha(matr) + +c = colorbar('Limits',[0 0.05]); +caxis([0 0.05]) +ax = gca; +ax.FontSize = 12; +ax.FontName = 'Arial'; +ax.XColor = [0 0 0]; +ax.YColor = [0 0 0]; +ax.XLabel.Color = [0 0 0]; +ax.YLabel.Color = [0 0 0]; +xtickangle(45) +xticks([1:6]) +yticks([1:6]) + +xticklabels({'M2','PPC','ACC', 'Thal', 'V2','OFC'}) +yticklabels({'M2','PPC','ACC', 'Thal', 'V2','OFC'}) + + +%% decoding vs cell no +for r = 1:6%length(area) + rmse_area{r} = []; lengthCells = [];cellsSpaced = []; rand_pos{r}= []; rmse_rand{r}= [];rmse_rand_CI{r}= []; + cellNo{r} = []; + for f = 1:length(file{r}) + +% file = files{r}{f}; + load(strcat('/Users/annachristinagarvert/Desktop/Michele/analysis/decoding/decodingOnlyPC/all_place_cells/',area{r},'/',file,'/d_data_varN.mat')) +% +% if r == 7 +% load(fullfile('/Users/annachristinagarvert/Desktop/Michele/sData/',area{r},'/',file,'/sData.mat')) +% else + load(fullfile('/Users/annachristinagarvert/Dropbox (UIO Physiology Dropbox)/Lab Data/Michele Gianatti/Data&Analysis/sData_SNR/',area{r},'/',file,'/final_sData.mat')) + +% end + + noTrials = length(trials.trialStart); + rmse = []; meanRMSE = []; mean_rand = []; rand_rmse = []; + for i = 1:length(d_data_varN) + + cellNo{r}{f}(i) = length(d_data_varN(i).iter{1, d_data_varN(i).medianTestFold}.cellIdx); + + +% for l = 1:5 + predPos = d_data_varN(i).iter{1, d_data_varN(i).medianTestFold}.predPos_test; + realPos = d_data_varN(i).iter{1, d_data_varN(i).medianTestFold}.realPos_test; + + maxval =154.3735;% max([max(predPos(:)) max(realPos(:))]); + predError = zeros(length(realPos),1); + for k = 1:length(realPos) + predError(k) = abs(realPos(k)-predPos(k)); + if predError(k) > maxval/2 && realPos(k) > predPos(k) + predError(k) = abs(maxval-realPos(k)+predPos(k)); + elseif predError(k) > maxval/2 && realPos(k) < predPos(k) + predError(k) = abs(maxval-predPos(k)+realPos(k)); + end + end + + meanRMSE(i) = nanmean(predError); + + for m = 1:5000 + + randPos = circshift(predPos, randperm(length(predPos),1)); +% if m == 10 +% figure() +% % plot(realPos); hold on +% plot(randPos) +% hold on +% plot(predPos) +% end + maxval =154.3735;% + randError = zeros(length(realPos),1); + for k = 1:length(realPos) + randError(k) = abs(realPos(k)-randPos(k)); + if randError(k) > maxval/2 && realPos(k) > randPos(k) + randError(k) = abs(maxval-realPos(k)+randPos(k)); + elseif randError(k) > maxval/2 && realPos(k) < randPos(k) + randError(k) = abs(maxval-randPos(k)+realPos(k)); + end + end + + rand_rmse(m,i) = nanmean(randError); + + if i == length(d_data_varN) + classes = unique(realPos); + realPosBinned = []; + classes = 2.6165:5.233:157.05; + for t = 1:length(realPos) + [~,ind] = min(abs(classes-realPos(t))); + realPosBinned(t) = classes(ind); + end + + + for k= 1:length(classes) + idx = find(realPosBinned == classes(k)); + rand_pos{r}(m,k) = nanmean(randError(idx)); + end + + end + + + + + end +% + +% % classes = unique(realPos); +% % if length(classes) > 30 +% % classes = [classes 149.1405+5.233]; +% % end +% +% classes = floor(2.6165:5.233:157.05); +% for k= 1:length(classes) +% idx = find(floor(realPos) == classes(k)); +% % mouseTemp(l,k) = nanmean(predError(idx)); +% mouse{r}(f,m,k) = nanmean(randError(idx)); +% end +% +% rand_rmse{r}(f,m) =nanmean(randError); +% +% end + end +% end + + mean_rand = nanmean(rand_rmse); + tempCI = quantile(rand_rmse,[0.05 0.95]); mean_rand_CI = tempCI(1,:); + if f > 1 + if size(rmse_area{r},2) < length(meanRMSE) + rmse_rand{r} = [rmse_rand{r} NaN(size(rmse_area{r},1),length(meanRMSE)-size(rmse_area{r},2))]; + rmse_rand_CI{r} = [rmse_rand_CI{r} NaN(size(rmse_area{r},1),length(meanRMSE)-size(rmse_area{r},2))]; + rmse_area{r} = [rmse_area{r} NaN(size(rmse_area{r},1),length(meanRMSE)-size(rmse_area{r},2))]; + end + if length(meanRMSE) < size(rmse_area{r},2) + meanRMSE = [meanRMSE NaN(1,size(rmse_area{r},2)-length(meanRMSE))]; + mean_rand = [mean_rand NaN(1,size(rmse_area{r},2)-length(mean_rand))]; + mean_rand_CI = [mean_rand_CI NaN(1,size(rmse_area{r},2)-length(mean_rand_CI))]; + end + end + + rmse_area{r}(f,:) = meanRMSE; + rmse_rand{r}(f,:) = mean_rand; + rmse_rand_CI{r}(f,:) = mean_rand_CI; + + cellsSpaced(f) = 10*round(cellNo{r}{f}(2)/10); + lengthCells(f) = length(cellNo{r}{f}); + end + + + if ~isempty(find(cellsSpaced == 20)) && ~isempty(find(cellsSpaced == 30)) + cellsSpacedShort = find(cellsSpaced == 20); + lengthCells(cellsSpacedShort) = 1; + for f = 1:length(cellsSpacedShort) + rmse_area{r}(cellsSpacedShort(f),2:2:end) = NaN; + end + end + [maxCells maxCellIdx]= max(lengthCells); + cellNoSum = cellNo{r}{maxCellIdx}; +% for f = 1:size(rmse_area{r},1) +% scatter(cellNo, rmse_area{r}(f,:) ,'k', 'filled') +% hold on +% end + noNan = []; + for m = 1:size(rmse_area{r},2) + noNan(m) = length(find(isnan(rmse_area{r}(:,m)))); + end + + idx = find(noNan == size(rmse_area{r},1)-1); + if ~isempty(idx) + plot(gca,cellNoSum(1:idx(1)-1),nanmean(rmse_area{r}(:,1:idx(1)-1),1), 'Color',col(r,:), 'LineWidth', 1.5) + hold(gca,'on') + plot(gca,cellNoSum(idx(1)-1:end),nanmean(rmse_area{r}(:,idx(1)-1:end),1), 'Color',col(r,:), 'LineWidth', 1.5, 'LineStyle','--') + xval = cellNoSum(1:idx(1)-1); yval = rmse_area{r}(:,1:idx(1)-1); +% errorbar(axes(4),xval, nanmean(yval), nanstd(yval)./sqrt(size(rmse_area{r},1)-noNan(1:idx(1)-1)), 'LineWidth' ,1.5, 'Color',col(r,:)) + patch(gca,[xval fliplr(xval)],[nanmean(yval)+nanstd(yval)./sqrt(size(rmse_area{r},1)-noNan(1:idx(1)-1)) fliplr(nanmean(yval)-nanstd(yval)./sqrt(size(rmse_area{r},1)-noNan(1:idx(1)-1)))],col(r,:),'linestyle','none','FaceAlpha', 0.2); + else + plot(gca, cellNoSum,nanmean(rmse_area{r}), 'Color',col(r,:), 'LineWidth', 1.5) + hold(gca,'on') + xval = cellNoSum; yval = rmse_area{r}; +% errorbar(axes(4),xval, nanmean(yval), nanstd(yval)./sqrt(size(rmse_area{r},1)), 'LineWidth' ,1.5, 'Color',col(r,:)) + patch(gca,[xval fliplr(xval)],[nanmean(yval)+nanstd(yval)./sqrt(size(rmse_area{r},1)) fliplr(nanmean(yval)-nanstd(yval)./sqrt(size(rmse_area{r},1)))],col(r,:),'linestyle','none','FaceAlpha', 0.2); + + end + +% errorbar(cellNo(1:idx(1)),nanmean( rmse_area{r}(:,1:idx(1)),1), nanstd( rmse_area{r}(:,1:idx(1)),1), 'k', 'LineStyle', 'none', 'LineWidth', 1.5) +% yline(nanmean(rand_rmse(:)), 'k', 'LineWidth', 1.5) + +% xlimval = [ 0 max(cellNo)+20 max(cellNo)+20 0]; +% ylimval = [nanmean(rand_rmse(:))-2*nanstd(rand_rmse(:)) nanmean(rand_rmse(:))-2*nanstd(rand_rmse(:)) nanmean(rand_rmse(:)) nanmean(rand_rmse(:))]; +% patch(xlimval, ylimval,[0.5 0.5 0.5],'LineStyle','none','FaceAlpha', 0.2); +% xlim([ 0 max(cellNo)+20]) + % yline(nanmean(rand_rmse(:))-2*nanstd(rand_rmse(:))) + + +% ylim([5 40]) +end + +xlabel(gca,'# axons') +ylabel(gca,'Decoding error (cm)') +set(gca, 'FontName', 'Arial', 'FontSize', 12) +box(gca, 'off') +ylim(gca,[5 45]) + + +figure() +dir = [2 4 3 1 5]; +for r = 2:6 + for l = 1:size(rmse_area{r},1) + scatter(dir(r-1),rmse_area{r}(l,14),70, [0.5 0.5 0.5]); + hold on + end + scatter(dir(r-1),nanmean(rmse_area{r}(:,14)),70, [0 0 0],'filled'); + errorbar(dir(r-1),nanmean(rmse_area{r}(:,14)),nanstd(rmse_area{r}(:,14)),'k','LineWidth',1.5); +end +xlim([0 6]) +ylim([5 55]) +xticks([1:5]) +xticklabels({'PPC','M2','Thal', 'ACC', 'V2'}) + +figure() +randall = []; +for r = 1:6 + randall = [randall ;rand_pos{r}]; + hold on +end +plot(nanmean(randall)) +tempCI = quantile(randall,[0.05 0.95]); +hold on +plot(tempCI(1,:)) + +rmse_rand_all = []; rmse_rand_CI_all = []; +for i = 1:6 + % figure() + for f = 1:length(cellNo{i}) + + v = rmse_rand_all; t = rmse_rand_CI_all; + if ~isempty(find(round(cellNo{i}{f}./10).*10 == 20)) + j = rmse_rand{i}(f,1:2:end); + k = rmse_rand_CI{i}(f,1:2:end); + else + j = rmse_rand{i}(f,:); + k = rmse_rand_CI{i}(f,:); + end + + n =max(size(v,2),numel(j)); + v(:,end+1:n)=nan; + j(end+1:n)=nan; + rmse_rand_all = [v; j]; + + n =max(size(t,2),numel(k)); + t(:,end+1:n)=nan; + k(end+1:n)=nan; + rmse_rand_CI_all = [t; k]; + end + % plot(cellNo{i}{f},rmse_rand{i}(f,~isnan(rmse_rand{i}(f,:)))) + % hold on +end +% hold on +% +% for i = 1:7 +% cellNoMax = length +% cell2mat(cellNo) +% +% rmse_rand{1} = rmse_rand{1}(:,1:2:end); +% for i = 1:7 +% for f = 1:length(cellNo{i}) +% if find(cellNo{i}{f} +% +% end + +F = 10:20:730; +amean = nanmean(rmse_rand_all(:,1:end-3)); +astd = nanmean(rmse_rand_CI_all(:,1:end-3)); +yline(nanmean(amean), 'LineWidth', 1.2, 'Color', 'k') +yline(nanmean(astd), 'LineWidth', 1.2,'LineStyle','--', 'Color', 'k') + +hold(axes(4), 'on') +patch(axes(4),[F(~isnan(amean)) fliplr(F(~isnan(amean)))],[amean(~isnan(amean)) fliplr(astd)],'k','linestyle','none','FaceAlpha', 0.2); + + +% meanAll_ax=axes(2); hold on +% for j = 1:3%length(mouse) +% % fig = figure(); +% % hax=axes; hold on +% +% if length(mouse{j}) > 1 +% xLength = []; x = []; y = []; yAll = []; +% +% for i = 1: length(mouse{j}) +% if mouseIdx{j}{i} == 1 +% uiopen(strcat(mouse{j}{i}, '/decoding_varN/rmse_varN.fig'),1); +% else +% uiopen(strcat(mouse{j}{i}, '/rmse_varN.fig'),1); +% end +% +% dataObjs = findobj(gca,'-property','YData'); +% % plot(hax,dataObjs(4).XData ,dataObjs(4).YData, 'LineWidth',2, 'Color', 'k');%col{j}) +% x{i} = dataObjs(4).XData; +% y{i} = dataObjs(4).YData; +% +% % if j ~= 4 || length(x{i})>1 +% if x{i}(2) == 20 +% y{i} = y{i}(1:2:end); +% x{i} = x{i}(1:2:end); +% end +% % end +% xLength = [xLength, length(x{i})]; +% end +% +% set(gca, 'FontName', 'Gotham', 'FontSize', 18) +% xlabel('number of axons') +% ylabel('test mean rmse [cm]') +% box off +% +% [maxLength maxI] = max(xLength); +% for i = 1: length(mouse{j}); yAll(i,:)= [y{i}, NaN(1,maxLength-xLength(i))]; end +% % patch(meanAll_ax,[x{maxI} fliplr(x{maxI})],[nanmean(yAll)+nanstd(yAll)./length(mouse{j}) fliplr(nanmean(yAll)-nanstd(yAll)./length(mouse{j}))],col(j,:),'linestyle','none','FaceAlpha', 0.4); +% plot(axes(4),x{maxI},nanmean(yAll), 'LineWidth' ,1.5, 'Color', col(j,:)) +% hold(axes(4),'on') +% errorbar(axes(4),x{maxI}, nanmean(yAll), nanstd(yAll), 'LineWidth' ,1.5, 'Color',col(j,:)) +% set(axes(4), 'FontName', 'Arial', 'FontSize', 8) +% xlabel(axes(4),'# of axons') +% ylabel(axes(4),'test mean rmse [cm]') +% box(axes(4),'off') +% else +% uiopen(strcat( mouse{j}{1}, '/decoding_varN/rmse_varN.fig'),1); +% dataObjs = findobj(gca,'-property','YData'); +% x{i} = dataObjs(4).XData; +% y{i} = dataObjs(4).YData; +% if x{i}(2) == 20 +% y{i} = y{i}(1:2:end); +% x{i} = x{i}(1:2:end); +% end +% % plot(hax,x{i},y{i}, 'LineWidth',2, 'Color', col(j,:)) +% +% box on +% set(gca, 'FontName', 'Gotham', 'FontSize', 22) +% xlabel('# cells') +% ylabel('test rmse [cm]') +% plot(axes(4),dataObjs(4).XData,dataObjs(4).YData, 'LineWidth' ,2, 'Color', col(j,:)) +% +% end +% end + + +function [corrected_p, h]=bonf_holm(pvalues,alpha) +if nargin<1 + error('You need to provide a vector or matrix of p-values.'); +else + if ~isempty(find(pvalues<0,1)), + error('Some p-values are less than 0.'); + elseif ~isempty(find(pvalues>1,1)), + fprintf('WARNING: Some uncorrected p-values are greater than 1.\n'); + end +end +if nargin<2 + alpha=.05; +elseif alpha<=0 + error('Alpha must be greater than 0.'); +elseif alpha>=1 + error('Alpha must be less than 1.'); +end +s=size(pvalues); +if isvector(pvalues) + if size(pvalues,1)>1 + pvalues=pvalues'; + end + [sorted_p sort_ids]=sort(pvalues); +else + [sorted_p sort_ids]=sort(reshape(pvalues,1,prod(s))); +end +[dummy, unsort_ids]=sort(sort_ids); %indices to return sorted_p to pvalues order +m=length(sorted_p); %number of tests +mult_fac=m:-1:1; +cor_p_sorted=sorted_p.*mult_fac; +cor_p_sorted(2:m)=max([cor_p_sorted(1:m-1); cor_p_sorted(2:m)]); %Bonferroni-Holm adjusted p-value +corrected_p=cor_p_sorted(unsort_ids); +corrected_p=reshape(corrected_p,s); +h=corrected_p 0 + + predPos = mData(r,l).d_data_onlyPC(end).iter{1, mData(r,l).d_data_onlyPC(end).medianTestFold}.predPos_test; + realPos = mData(r,l).d_data_onlyPC(end).iter{1, mData(r,l).d_data_onlyPC(end).medianTestFold}.realPos_test; + + maxval =154.3735;% max([max(predPos(:)) max(realPos(:))]); + predError = zeros(length(realPos),1); + for k = 1:length(realPos) + predError(k) = abs(realPos(k)-predPos(k)); + if predError(k) > maxval/2 && realPos(k) > predPos(k) + predError(k) = abs(maxval-realPos(k)+predPos(k)); + elseif predError(k) > maxval/2 && realPos(k) < predPos(k) + predError(k) = abs(maxval-predPos(k)+realPos(k)); + end + end + + rmse(r,f) = nanmean(predError); + f = f+1; + end + end + +end + +rmse(rmse==0) = NaN; + +%% decoding errors for each area +fig_3_e = figure(); + +sorting = [6 1 3 4 2 5]; +for r = 1:6 + for l = 1:size(rmse,2) + scatter(sorting(r),rmse(r,l),70, [0.5 0.5 0.5]); + hold on + end + scatter(sorting(r),nanmean(rmse(r,:)),70, [0 0 0],'filled'); + errorbar(sorting(r),nanmean(rmse(r,:)),nanstd(rmse(r,:))/sqrt(sum(~isnan(rmse(r,:)))),'k','LineWidth',1.5); +end + +xlim([0 7]) +ylim([5 55]) +xticks([1:6]) +ax = gca; +ax.FontSize = 12; +ax.FontName = 'Arial'; +ax.XColor = [0 0 0]; +ax.YColor = [0 0 0]; +ax.XLabel.Color = [0 0 0]; +ax.YLabel.Color = [0 0 0]; +xticklabels({'M2','PPC','ACC', 'Thal', 'V2M','OFC'}) + + +sorting_inverse =[2 5 3 4 6 1]; + +for i = 1:6 + for j = 1:6 + pval(i,j) = ranksum(rmse(sorting_inverse(i),:),rmse(sorting_inverse(j),:)); + end +end + +matr =triu(ones(6,6)); +matr(eye(6,6)==1) = NaN; +matr(matr== 0)= NaN; +pval(isnan(matr)) = NaN; +[corrected_p, h]=helper.bonf_holm(pval); + + +%% plot p values +fig_3_f = figure(); +imagesc(triu(corrected_p)) +hold on +alpha(matr) + +c = colorbar('Limits',[0 0.05]); +caxis([0 0.05]) +ax = gca; +ax.FontSize = 12; +ax.FontName = 'Arial'; +ax.XColor = [0 0 0]; +ax.YColor = [0 0 0]; +ax.XLabel.Color = [0 0 0]; +ax.YLabel.Color = [0 0 0]; +xtickangle(45) +xticks([1:6]) +yticks([1:6]) + +xticklabels({'M2','PPC','ACC', 'Thal', 'V2','OFC'}) +yticklabels({'M2','PPC','ACC', 'Thal', 'V2','OFC'}) + +% \ No newline at end of file diff --git a/Figures/figure_3/fig_3_a.m b/Figures/figure_3/fig_3_a.m new file mode 100644 index 0000000..f3817dd --- /dev/null +++ b/Figures/figure_3/fig_3_a.m @@ -0,0 +1,34 @@ + + +% OFC +area = 1; +session = 8; +plot.plotRealVsDecodedPosition(mData,area,session) + +% M2 +area = 2; +session = 8; +plot.plotRealVsDecodedPosition(mData,area,session) + +% ACC +area = 3; +session = 10; +plot.plotRealVsDecodedPosition(mData,area,session) + +% Thalamus +area = 4; +session = 10; +plot.plotRealVsDecodedPosition(mData,area,session) + +% PPC +area = 5; +session = 3; +plot.plotRealVsDecodedPosition(mData,area,session) + + +% V2 +area = 6; +session = 7; +plot.plotRealVsDecodedPosition(mData,area,session) + + diff --git a/Figures/figure_3/fig_3_d_g.m b/Figures/figure_3/fig_3_d_g.m new file mode 100644 index 0000000..3eeed11 --- /dev/null +++ b/Figures/figure_3/fig_3_d_g.m @@ -0,0 +1,28 @@ +col = [122 172 210 + 80 106 139 + 217 85 88 + 129 89 162 + 180 151 94 + 179 205 142 + 0 0 0 + 255 50 80]./255; + +fig_d = figure(); + +for a = 1:length(data)-2 + area = a; + plot.decodingVsNoCells + + % plot results + plot(gca,cellNoSum,nanmean(rmse_all,1), 'Color',col(area,:), 'LineWidth', 1.5) + yval = rmse_all; + patch(gca,[cellNoSum fliplr(cellNoSum)],[nanmean(yval)+nanstd(yval)./sqrt(size(rmse_all,1)) fliplr(nanmean(yval)-nanstd(yval)./sqrt(size(rmse_all,1)))],col(area,:),'linestyle','none','FaceAlpha', 0.2); + hold on +end + + +ylabel('Decoding error (cm)') +xlabel('# Cells') +set(gca,'FontName','Arial','FontSize',12) +box off +ylim([0 55]) \ No newline at end of file diff --git a/Figures/figure_3/fig_3_h.m b/Figures/figure_3/fig_3_h.m new file mode 100644 index 0000000..ca7e731 --- /dev/null +++ b/Figures/figure_3/fig_3_h.m @@ -0,0 +1,85 @@ +% the skaggs values were calculated with the function % calculate skaggs + +SI = cell(length(data),1); +SIPC = cell(length(data),1); +SI_notPC{l} = cell(length(data),1); + +for l = 1:length(data) + SIPC{l} = []; + + for i = 1: length(data(l).sessionIDs) + + SI_session = mData(l,i).skaggs(:); + + % find all place tuned cells + pcNo = []; + for k = 1:length(mData(l,i).LNP_results) + if ~isempty(find(mData(l,i).LNP_results(k).selected_model == 1)) && isempty(find(mData(l,i).LNP_results(k).selected_model == 5)) + pcNo = [pcNo k]; + end + end + + SIPC{l} = [SIPC{l}; real(SI_session(pcNo))]; + meanSIPC{l}(i) = nanmean(real(SI_session(pcNo))); + + end +end + + +figure() +subplot(1,2,1) +violin({SIPC{2},SIPC{5},SIPC{4},... + SIPC{3},SIPC{6},SIPC{1}},'facecolor',[0.5 0.5 0.5]) +xticklabels({ 'M2', 'PPC', 'Thal', 'ACC', 'V2', 'OFC'}) +xticks(gca,[1:6]) +ylabel('Skaggs Spatial Information (Bits/AP)') +ax = gca; +ylim([0 12]) +ax.FontSize = 12; +ax.FontName = 'Arial'; +ax.XColor = [0 0 0]; +ax.YColor = [0 0 0]; +ax.XLabel.Color = [0 0 0]; +ax.YLabel.Color = [0 0 0]; +box off + + +sortVal = [2 5 4 3 6 1]; + +% calculate all pvals with ranksum +for i = 1:6 + for j = 1:6 + pval(i,j) = ranksum(meanSIPC{sortVal(i)}',meanSIPC{sortVal(j)}'); + end +end + +matr =triu(ones(6,6)); +matr(eye(6,6)==1) = NaN; +matr(matr== 0)= NaN; +pval(isnan(matr)) = NaN; +[corrected_p, ~]=helper.bonf_holm(pval); + +subplot(1,2,2) +imagesc(triu(corrected_p)) +hold on +alpha(matr) + +xticks([1:6]) +yticks([1:6]) +c = colorbar('Limits',[0 0.05]); +caxis([0 0.05]) + +xticklabels({'M2', 'PPC', 'Thal', 'ACC', 'V2', 'OFC'}) +yticklabels({'M2', 'PPC', 'Thal', 'ACC', 'V2', 'OFC'}) + +ax = gca; +ax.FontSize = 12; +ax.FontName = 'Arial'; +ax.XColor = [0 0 0]; +ax.YColor = [0 0 0]; +ax.XLabel.Color = [0 0 0]; +ax.YLabel.Color = [0 0 0]; +xtickangle(45) +title('Place tuned axons/cells') + + diff --git a/Figures/load_LNP_results.m b/Figures/load_LNP_results.m new file mode 100644 index 0000000..d25c3d4 --- /dev/null +++ b/Figures/load_LNP_results.m @@ -0,0 +1,11 @@ + + +LNP_dirct = '/Users/annachristinagarvert/UIO Physiology Dropbox Dropbox/Lab Data/Michele Gianatti/Data&Analysis/Paper/LNP'; + +for f = 1:length(data) + for i = 1:length(data(f).sessionIDs) + load(fullfile(LNP_dirct,data(f).area,strcat(data(f).sessionIDs{i},'.mat'))) + mData(f,i).LNP_results = results; + end +end + diff --git a/Figures/load_all_data.m b/Figures/load_all_data.m new file mode 100644 index 0000000..bf918e3 --- /dev/null +++ b/Figures/load_all_data.m @@ -0,0 +1,129 @@ +%% Areas to include +data(1).area = 'OFC'; +data(2).area = 'M2'; +data(3).area = 'ACC'; +data(4).area = 'Thal'; +data(5).area = 'PPC'; +data(6).area = 'V2'; +data(7).area = 'RSC'; +data(8).area = 'V1'; + + +%% Recordings to include for area 1-8 +data(1).sessionIDs = {'m5143-20210710-02',... + 'm5142-20210709-02',... + 'm5142-20210709-01',... + 'm5112-20201202-02',... + 'm5112-20201125-02',... + 'm5112-20201108-01',... + 'm5111-20201202-01',... + 'm5111-20201125-02',... + 'm5111-20201125-01',... + 'm5110-20201125-03',... + 'm5110-20201125-01',... + 'm5110-20201202-01',... + 'm5110-20201108-01'}; + +data(2).sessionIDs = {'m5094-20200520-01',... + 'm5094-20200522-01',... + 'm5094-20200522-03',... + 'm5095-20200522-01',... + 'm5095-20200520-01',... + 'm5122-20210304-01',... + 'm5122-20210305-01',... + 'm5128-20210324-01'}; + +data(3).sessionIDs = {'m5125-20210325-01',... + 'm5125-20210324-01',... + 'm5125-20210323-01',... + 'm5124-20210324-02',... + 'm5124-20210323-02',... + 'm5085-20200314-01',... + 'm5085-20200313-02',... + 'm5085-20200313-01',... + 'm5082-20200313-02',... + 'm5082-20200312-02',... + 'm5082-20200312-01'}; + +data(4).sessionIDs = {'m5074-20191202-01',... + 'm5075-20191203-01',... + 'm5074-20191205-01',... + 'm5074-20191219-01',... + 'm5075-20191205-01',... + 'm5075-20191205-02',... + 'm5075-20191219-01',... + 'm5088-20200430-02',... + 'm5088-20200430-03',... + 'm5088-20200504-01',... + 'm5088-20200504-02',... + 'm5088-20200504-03',... + 'm5088-20200508-01',... + 'm5088-20200508-03',... + 'm5140-20210709-01',... + 'm5140-20210710-01',... + 'm5140-20210710-02',... + 'm5141-20210709-02',... + 'm5141-20210710-01'}; + + +data(5).sessionIDs = {'m5119-20210204-02',... + 'm5121-20210204-03',... + 'm5121-20210205-01',... + 'm5121-20210205-02',... + 'm5119-20210204-01',... + 'm5150-20211207-02',... + 'm5150-20211208-01',... + 'm5153-20211208-01',... + 'm5153-20211208-02'}; + +%% V2 +data(6).sessionIDs = {'m5137-20210610-01',... + 'm5137-20210610-02',... + 'm5138-20210610-01',... + 'm5138-20210610-02',... + 'm5138-20210611-01',... + 'm5139-20210611-01',... + 'm5139-20210611-02',... + 'm5139-20210612-01',... + 'm5146-20211021-01',... + 'm5146-20211024-02',... + 'm5146-20211022-01',... + 'm5147-20211022-01',... + 'm5147-20211024-01'}; + +%% RSC +data(7).sessionIDs = {'m5115-20201203-02',... + 'm5115-20201203-01',... + 'm5115-20201202-01',... + 'm5114-20201203-01',... + 'm5114-20201202-02',... + 'm5114-20201202-01',... + 'm5108-20201120-03',... + 'm5108-20201120-02',... + 'm5108-20201120-01'}; + + +data(8).sessionIDs = {'m5091-20200504-03',... + 'm5091-20200508-03',... + 'm5091-20200508-06',... + 'm5092-20200520-01',... + 'm5092-20200603-01',... + 'm5092-20200603-03',.... + 'm5101-20200714-01',... + 'm5102-20200713-01',... + 'm5102-20200713-02',... + 'm5102-20200714-02'}; + + + +load_sData_dirct = '/Users/annachristinagarvert/UIO Physiology Dropbox Dropbox/Lab Data/Michele Gianatti/Data&Analysis/Paper/sData/'; +analysis_dirct = '/Users/annachristinagarvert/UIO Physiology Dropbox Dropbox/Lab Data/Michele Gianatti/Data&Analysis/Paper/'; + +mData = struct(); +% Load all sessions +for i = 1:length(data) + for f = 1:length(data(i).sessionIDs) + mData(i,f).sData= load(fullfile(load_sData_dirct,data(i).area,data(i).sessionIDs{f}, 'final_sData.mat')); + end +end + diff --git a/Figures/load_decoding_results.m b/Figures/load_decoding_results.m new file mode 100644 index 0000000..15d4ce7 --- /dev/null +++ b/Figures/load_decoding_results.m @@ -0,0 +1,20 @@ + +decoding_onlyPC_dirct = '/Users/annachristinagarvert/UIO Physiology Dropbox Dropbox/Lab Data/Michele Gianatti/Data&Analysis/Paper/decoding/decodingOnlyPC/'; + +for f = 1:length(data) + for i = 1:length(data(f).sessionIDs) + load(fullfile(decoding_onlyPC_dirct,data(f).area,data(f).sessionIDs{i},'d_data_varN.mat')) + mData(f,i).d_data_onlyPC = d_data_varN; + end +end + +decoding_dirct = '/Users/annachristinagarvert/UIO Physiology Dropbox Dropbox/Lab Data/Michele Gianatti/Data&Analysis/Paper/decoding/decodingVsNoCells/'; + +for f = 1:length(data) + for i = 1:length(data(f).sessionIDs) + load(fullfile(decoding_dirct,data(f).area,data(f).sessionIDs{i},'d_data_varN.mat')) + mData(f,i).d_data = d_data_varN; + end +end + + diff --git a/Figures/load_skaggs_results.m b/Figures/load_skaggs_results.m new file mode 100644 index 0000000..f17f37f --- /dev/null +++ b/Figures/load_skaggs_results.m @@ -0,0 +1,10 @@ + +decoding_skaggs = '/Users/annachristinagarvert/UIO Physiology Dropbox Dropbox/Lab Data/Michele Gianatti/Data&Analysis/Paper/Skaggs/'; + +for f = 1:length(data) + for i = 1:length(data(f).sessionIDs) + load(fullfile(decoding_skaggs,data(f).area,data(f).sessionIDs{i},'SI_session.mat')) + mData(f,i).skaggs = SI_session; + end +end + diff --git a/Figures/sData_template.m b/Figures/sData_template.m new file mode 100644 index 0000000..65042cd --- /dev/null +++ b/Figures/sData_template.m @@ -0,0 +1,172 @@ +%% The following nomenclature must be used for sData +% Updated 02-10-2018 + + +%% sData +sData = struct(); % Initialize the struct + +sData.mouseInfo = struct(); % Struct containing all information about the mouse +sData.sessionInfo = struct(); % Struct containing information about the session +sData.imdata = struct(); % Imaging data +sData.daqdata = struct(); % Signals obtained from the control computer DAQ +sData.ephysdata = struct(); % Ephys signals +sData.behavior = struct(); % Behavioral variables such as for example lick rate +sData.trials = struct(); % Useful if the experiment contains trials. +sData.stats = struct(); % Useful to add analysis variables and statistics + + +%% MOUSEINFO +sData.mouseInfo.name % char (REQUIRED) Mouse number according to lab convention +sData.mouseInfo.dateOfBirth % char (REQUIRED) yyyy.mm.dd +sData.mouseInfo.strain % char (REQUIRED) Strain of animal used. Be specific. +sData.mouseInfo.sex % char (REQUIRED) Male or female +sData.mouseInfo.surgeryDate % char (REQUIRED) yyyy.mm.dd +sData.mouseInfo.windowCoordinates % double (REQUIRED?) [x, y] (mm) distance from bregma to center of window. +sData.mouseInfo.lightCycle % char (Optional) Light is on during the following time, ex: '10 am to 10 pm'. + +sData.mouseInfo.SLCageNumber % double (Optional) Science linker cage number +sData.mouseInfo.SLMouseNumber % double (Optional) Science Linker animal number +sData.mouseInfo.surgeryProtocol % char (Optional) should refer to a documented protocol +sData.mouseInfo.injectedVirus % char (Optional) Relevant for injections +sData.mouseInfo.hemisphere % char (Optional) relevant for injections/recordings +sData.mouseInfo.injectionCoordinates % double (Optional) nInj x 2 (um) from bregma + + +%% SESSIONINFO +sData.sessionInfo.sessionID % char (REQUIRED) A unique ID for the session. It takes the form : mouseName-YYYYMMDD-SS (where SS is the session number). +sData.sessionInfo.date % char (REQUIRED) yyyy.mm.dd +sData.sessionInfo.sessionNumber % double (REQUIRED) The session number part of your sessionID. +sData.sessionInfo.sessionStartTime % char (REQUIRED) hh:mm:ss +sData.sessionInfo.sessionStopTime % char (REQUIRED) hh:mm:ss +sData.sessionInfo.recordedData % cell (REQUIRED): Either {'2P','LFP','Patch','PupilVideo','SecurityCamera'}. These are only peripheral recordings for your DAQ, i.e. you do not need to add running wheel etc. + +sData.sessionInfo.mouseWeight % double (REQUIRED, for water deprivation experiments) +sData.sessionInfo.mouseOriginalWeight % double (REQUIRED, for water deprivation experiments) +sData.sessionInfo.wasSessionAborted % logic (Optional) Was the session aborted before you wanted it to? +sData.sessionInfo.recDayNumber % double (Optional) 1,2,3,...,N +sData.sessionInfo.experimentName % char (Optional) name of experiment +sData.sessionInfo.protocol % char (Optional) Tells where a documented protocol of your experiments are located. Fex Google docs etc. + + +%% DAQDATA +sData.daqdata.lickSignal % double Lick signal +sData.daqdata.wheelRotaryEncoderSignal % double Rotary encoder signal used for running wheels +sData.daqdata.wheelDiode % double Photodiode signal used to signify absolute position of wheel +sData.daqdata.waterValve % double Water valve signal showing when the water valve is open +sData.daqdata.frameSignal % double Frame signal recorded from the microscope +sData.daqdata.optoSignal % double Optogenetics signal which is the voltage also sent to the laser diode +sData.daqdata.frameIndex % double Array of same lenght as frames required, where each sample n is the index number of other daqdata samples at the onset of frame n. + + +% DAQ Metadata +sData.daqdata.meta.fs % double (REQUIRED) Sampling frequency of the data acquisition system +sData.daqdata.meta.wheelRotaryEncoderTicks % double (REQUIRED, if rotary encoder wheel is used) Number of ticks on the rotary encoder wheel used + + +%% IMDATA +sData.imdata.roiSignals = struct(); % struct roiSignals is a struct array containing all signals as subfields. The channel color is specified as the array element number. +sData.imdata.roiSignals(2).ch % char (REQUIRED) Use roiSignals(x).ch to signify the color of the channel. This is different between OsloScope I and II. Specify channel number in () and add .ch as 'red' or 'green'. + +sData.imdata.roiSignals(2).roif % single ROI fluroescence +sData.imdata.roiSignals(2).npilf % single Neuropil fluorescence +sData.imdata.roiSignals(2).dff % single Delta F/F0 + +sData.imdata.roiSignals(2).deconv % single Deconvolved signal +sData.imdata.roiSignals(2).spikes % single Estimate spikes +sData.imdata.roiSignals(2).roifSubtractedNpil % single ROI fluorescence after subtracting neuropil +sData.imdata.roiSignals(2).dffSubtractedNpil % single Delta F/F0 after subtracting neuropil + +sData.imdata.time % single (Optional) If you want a time vector. This vector contains time stamps for each sample in roiSignals, starting at zero. + + +% Imaging Metadata (The following metadata can be import from the SciScan .ini file using the function getSciScanMetaData().) +sData.imdata.meta.microscope % char (REQUIRED) Name of microscope +sData.imdata.meta.xpixels % double (REQUIRED) Pixels in x dim +sData.imdata.meta.ypixels % double (REQUIRED) Pixels in y dim +sData.imdata.meta.dt % double (REQUIRED) Time difference between frames +sData.imdata.meta.fps % double (REQUIRED) Frames per second (average) +sData.imdata.meta.zoomFactor % double (REQUIRED) Zoom factor used for recording +sData.imdata.meta.zPosition % double (REQUIRED) Relative depth of imaging in Z dim. The experimenter zeros the relative depth to dura. +sData.imdata.meta.umPerPxX % double (REQUIRED) Micrometers per pixel in x dim +sData.imdata.meta.umPerPxY % double (REQUIRED) Micrometers per pixel in y dim +sData.imdata.meta.nChannels % double (REQUIRED) Number of recorded channels +sData.imdata.meta.channelNumbers % double (REQUIRED) Channels present, ex: [1,2]. +sData.imdata.meta.channelNames % cell (REQUIRED) Channel names, ex: {'Ch1','Ch2'}. +sData.imdata.meta.channelColor % cell (REQUIRED) Channel colors, ex: {'Red','Green'}. +sData.imdata.meta.nFrames % double (REQUIRED) Number of frames aquired +sData.imdata.meta.pockelSetting % double (REQUIRED) Pockel setting in SciScan in percent +sData.imdata.meta.pmtGain % double (REQUIRED) The gain set to PMT(s) during recording. (If two channels, [10, 20]) +sData.imdata.meta.piezoActive % logical (REQUIRED) true or false +sData.imdata.meta.piezoMode % char (REQUIRED) Description of the mode used on the piezo, either SAW or ZIG +sData.imdata.meta.piezoNumberOfPlanes % double (REQUIRED) Number of planes used by the piezo. +sData.imdata.meta.piezoImagingRateHz % double (REQUIRED) The volume frames per second rate (basically overall imaging rate divided by number of piezo planes). +sData.imdata.meta.piezoVolumeDepth % double (REQUIRED) The number of micrometer the volume imaging spans. + + +% The below has to be handwritten by experimenter +sData.imdata.meta.laserPower % double (REQUIRED) Laser power in mw +sData.imdata.meta.waveLength % double (REQUIRED) Wavelength used on the laser +sData.imdata.meta.fovCoordinates % double (REQUIRED) [x, y, theta] micrometer between center of FOV and center of window. Theta is the rotation in degree. + + +% roiArray +sData.imdata.roiArray % struct (REQUIRED) Struct array of same length as the number of ROIs. This is generated by roimanager and contains all metadata about the drawn ROIs, such as position in FOV, size etc. + + +%% EPHYSDATA +sData.ephysdata.lfp = []; + + +% Meta +sData.ephysdata.meta.samplingRate % double (REQUIRED) Sampling rate of ephys data +sData.ephysdata.meta.filterRange = [1000,100000]; % double (REQUIRED) Upper and lower +sData.ephysdata.meta.amplifier % char (REQUIRED) Name and type of amplifier +sData.ephysdata.meta.amplifierGain % double (REQUIRED) Gain on amplifier used +sData.ephysdata.meta.coordinateElectrodeTip % double (REQUIRED) [a.p, m.l., depth] +sData.ephysdata.meta.coordinateReferenceElectrodeTip % double (REQUIRED) [a.p, m.l., depth] + + +%% BEHAVIOR +% The naming of fields in behavior are as of yet up to each and every one +% to decide, but they should be perfectly explainable by reading the +% variable name. However, we advice you to use the following as +% abbreviations if needed: +% +% Ds +% All fields have a potential *Ds ending, which is ment for behavior +% downsampled to the frames acquired from the microscope or potentially +% some other signal, like LFP, tetrode recording, animal tracking etc. +% Binned +% Binned may refer to spatial or temporal binning. In the case of a running +% wheel, binning refers to spatial binning on the wheel with binSize being +% in cm. +% +% Examples follow: + +sData.behavior.wheelPos % double Absolute position on the wheel +sData.behavior.wheelPosDs % double Absolute position on the wheel (downsampled to the imaging FPS) +sData.behavior.wheelPosDsBinned % double Absolute position on the wheel (downsampled to the imaging FPS, and binned to some bin size stored in behavior.meta.binSize) + +sData.behavior.wheelLap % double Absolute lap number from the start of the recording. + +sData.behavior.runSpeed % double Running speed of the animal +sData.behavior.runSpeedDs % double Running speed of the animal (downsampled to the imaging FPS) + +% Meta +sData.behavior.meta.binSize % double Bin size when binning data. For instance, this can be the spatial bin size used for binning linear track data. +sData.behavior.meta.nWheelBins % double Number of bins a running wheel has been divided into (useful for place field recordings) + + +%% TRIALS +% The fields of trials are up to everyone to decide. But make sure that the +% variable names are logical and easy to understand. +sData.trials + + +%% STATS +% The field names used in stats are up to each individual to decide. But it +% is a conventient way of storing variables while you do analysis. +sData.stats + + + diff --git a/README.md b/README.md index b884b52..8a90334 100644 --- a/README.md +++ b/README.md @@ -1 +1,8 @@ -# RSC-projections \ No newline at end of file +# projections_to_RSC +Gianatti et. al,2023 + +Run "create_all_figures.m" in "Figures" to reproduce paper figures. +Functions require saved analysis. + +To replicate analysis use functions in "analysis". + diff --git a/analysis/cellClassification/calculateLLHIncrease.m b/analysis/cellClassification/calculateLLHIncrease.m new file mode 100644 index 0000000..df83532 --- /dev/null +++ b/analysis/cellClassification/calculateLLHIncrease.m @@ -0,0 +1,22 @@ +function sumMat = calculateLLHIncrease(signal,stateMatrix,param) + + % calculate log-likelihood increase from "mean firing rate model" + lambdaModel = exp(stateMatrix * param); + n = signal; + lambdaMean = nanmean(signal); + + LLHModel = exp(n.*log(lambdaModel) - lambdaModel - gammaln(n+1)); + log_LLHModel = nansum(-log(LLHModel))/sum(n); + + LLHMean = exp(n.*log(lambdaMean) - lambdaMean - gammaln(n+1)); + log_LLHMean = nansum(-log(LLHMean))/sum(n); + + logLLH_Increase = log(2)*(-log_LLHModel + log_LLHMean); + + sumMat = [logLLH_Increase log_LLHModel log_LLHMean]; + + +end + + + diff --git a/analysis/cellClassification/cellClassification.m b/analysis/cellClassification/cellClassification.m new file mode 100644 index 0000000..678a032 --- /dev/null +++ b/analysis/cellClassification/cellClassification.m @@ -0,0 +1,77 @@ + +%% Clear the workspace and load the data + +clear all; close all; clc +% load the data + +file = 'm5082-20200313-02'; +savePath = fullfile('/Users/annachristinagarvert/Desktop/',file); +if ~exist(savePath, 'dir'); mkdir(savePath); end + +sData = load(strcat('/Users/annachristinagarvert/Desktop/',file,'/final_sData.mat')); + +%% hyperparameters (these were optimised before) +N_P = 40; beta{1} = 0.01; % position +N_S = 30; beta{2} = 0.01; speedLimits = [0 60];% speed +N_A = 40; beta{3} = 0; accLimits = [-200 200];% acceleratioon +N_L = 20; beta{4} = 0; lickLimits = [0 10];% lick + +% other variables +noCov = 4; +numFold = 5; % attention this is also hard coded in the function fit_model +samplingRate = 10000; + +lickLengthThreshold = 15; +lickFreqThreshold = 15; +dt = 1/31; + +% prepare Signal +signal = double(sData.imdata.roiSignals(2).deconv); +numNeuron = size(deconvSignal,1); + +% prepare covariates +position = sData.behavior.wheelPosDs; +speed = sData.behavior.runSpeedDs; +acceleration = diff(runSpeed)./dt; +lickRate = extractLickFrequency(sData, samplingRate,lickLengthThreshold, lickFreqThreshold); + + +% make sure covariates and signal have the same length +minLength = min([length(position(2:end)),length(speed(2:end)),length(acceleration),length(lickRate(2:end)),size(signal(:,2:end),2)]); +signal = signal(:,2:minLength); + +position = position(2:minLength); +speed = speed(2:minLength); +acceleration = acceleration(1:minLength); +lickRate = lickRated(1:minLength); + +% create state matrix : +covGrid{1} = stateMatrixPosition(position, 0, max(position),N_P); % position +covGrid{2} = stateMatrix(speed, speedLimits(1), speedLimits(2),N_S); +covGrid{3} = stateMatrix(acceleration, accLimits(1), accLimits(2),N_A); +covGrid{4} = stateMatrix(lickRate, lickLimits(1), lickLimits(2),N_L); + +% calculate LLH +axon = cell(numNeuron,1); +selected_model = cell(numNeuron,1); + +for i = 1:numNeuron + for j = 1:noCov + [axon(i).testFit{j},axon(i).trainFit{j},axon(i).param{j}] = fit_model(covGrid{j},signal(i,:), beta{j}); + end + + axon(i).selected_model = compareModelPerformance(axon(i), numFold,noCov); + + [isCue, axon(i).cueData] = checkForCueTuning(sData,signal(i,:)); + + if isCue == 1;axon(i).selected_model = 5; end + +end + +LNP.axon = axonPos; save(fullfile(savePath,'axon'),'-struct', 'LNP'); + + + + + + \ No newline at end of file diff --git a/analysis/cellClassification/checkForCueTuning.m b/analysis/cellClassification/checkForCueTuning.m new file mode 100644 index 0000000..3a701d8 --- /dev/null +++ b/analysis/cellClassification/checkForCueTuning.m @@ -0,0 +1,63 @@ +function [isCue, cueData] = checkForCueTuning(sData, signal) +% Uses the signal (deconvolved) for a single cell and the sData struct +% and detemines wether the signal is cue tuned by first determiningthe place +% fields are determined and +% Then for cells that have more than one place field the distance between place +% fields is compared to the distance of the cues on the wheel + +% INPUT +% signal: Array of deconvolved signal for one cell with the length of number time bins + +% OUTPUT +% cueData: A struct containing all related cue field information +% .ifv: mean in-field value +% .field_width: width of place field as number of bins +% .field_peak_ind: index of peak of place field(s) +% .field_start: first index of place field(s) +% .field_end: last index of place field(s) +% .n_peaks: number of peaks in position activity map +% .inFieldIdxs: array containing all indices that are within any field of the position map +% .outOfFieldIdxs: array containing all indices that are considered out of a field +% .isCottonCell : boolean == 1 if the cell is a cotton cue cell +% (otherwise == 0) +% .isVelcroCell : boolean == 1 if the cell is a velcro cue cell (otherwise == 0) +% isCue : 1 if signal is cue tuned otherwise 0 + + +signalMatrix = splitInTrials(signal, sData); +smoothedMatrix = smoothdata(signalMatrix , 2,'gaussian',5); +meanSignal = nanmean(smoothedMatrix, 1); +cueData.fieldInfo = findPlaceFields(meanSignal,signalMatrix , sData); + +% if axon has more than one peak calculate peak distance +if cueData.fieldInfo.nPeaks >= 2 + peak = zeros(cueData.fieldInfo.nPeaks,1); + for m = 1:cueData.fieldInfo.nPeaks + peak(m) = cueData.fieldInfo.field_peak_ind(m); + end + + peakDiff = diff(sort(peak)); + + % position of cues: + % cue1 = [42 47 122 128]./2; + % cue2 = [63 68 100 106]./2; + % distance of peaks to be counted as cue tuned: + cottonDist = 16:22; % cue 2 : distance 19 +/- 3; + velcroDist = 37:43; % cue 1 : distance 80 +/-3; + + isCue = 0; + for m = 1:length(peakDiff) + if ~isempty(find(cottonDist == peakDiff(m), 1)) + cueData.isCottonCell = 1; + isCue = 1; + end + if ~isempty(find(velcroDist == peakDiff(m), 1)) + cueData.isVelcroCell = 1; + isCue = 1; + end + end +else + isCue = 0; +end + +end \ No newline at end of file diff --git a/analysis/cellClassification/compareModelPerformance.m b/analysis/cellClassification/compareModelPerformance.m new file mode 100644 index 0000000..10da1fa --- /dev/null +++ b/analysis/cellClassification/compareModelPerformance.m @@ -0,0 +1,46 @@ +function results = compareModelPerformance(results, numFold,noCov) + +testFit_mat = cell2mat(results.testFit); +LLH_values = reshape(testFit_mat(:,1:3:end),numFold,noCov); + +% find all models that performe above baseline +aboveBaseline = []; +for i = 1:noCov + pval_baseline = signrank(LLH_values(:,i),[],'tail','right'); + if pval_baseline <= 0.1 + aboveBaseline(end+1) = i; + end +end + +% check wether one model outperforms another if more than one model +% performe above baseline +sortOut = []; +if length(aboveBaseline) > 1 + for i = 1:length(aboveBaseline) + for j = 1:length(aboveBaseline) + [pval(i,j),~] = signrank(LLH_values(:,i),LLH_values(:,j),'tail','right'); + if pval(i,j) <= 0.1 % model i is significanlty better than model j + sortOut(end+1) = j; + end + end + end +end + +results.selectedModels = setdiff(aboveBaseline,sortOut); + +if ~isempty(results.selectedModels) + results.selectedModels = NaN; +end + +% save p-values for all models +if (sum(isnan(LLH_values(:,1))) == 5); results.P_pval = 1; else + results.P_pval = signrank(LLH_values(:,1),[],'tail','right'); end +if (sum(isnan(LLH_values(:,2))) == 5); results.S_pval = 1; else + results.S_pval = signrank(LLH_values(:,2),[],'tail','right'); end +if (sum(isnan(LLH_values(:,3))) == 5); results.A_pval = 1; else + results.A_pval = signrank(LLH_values(:,3),[],'tail','right');end +if (sum(isnan(LLH_values(:,4))) == 5); results.L_pval = 1; else + results.L_pval = signrank(LLH_values(:,4),[],'tail','right');end + + +end \ No newline at end of file diff --git a/analysis/cellClassification/findPlaceFields.m b/analysis/cellClassification/findPlaceFields.m new file mode 100644 index 0000000..317ccb0 --- /dev/null +++ b/analysis/cellClassification/findPlaceFields.m @@ -0,0 +1,228 @@ +function pfData = findPlaceFields(avgActivity,signalMatrix,sData) +% FINDPLACEFIELDS +% Uses the averaged position activity map (avgActivity) to extract information about +% potential place fields within the map. +% The function takes in an array of length corresponding to the number of +% bins on the wheel (ex 105), where each data point is the summed and +% position occupancy corrected average activity for all laps. A gaussian +% smoothing should also have been applied (like Mao et al (2017) does). +% The signal matrix activity in position bins x trials is further used to +% check if the place field fulfill the following criteria: +% --1 Criteria: the activityRatio of in field activity vs out of field activity is +% above 3 (ActivityRatioCriteria = 3) % This is the ratio used by Mao et al (2017) +% --2 Criteria: The percentage of laps where the activity peak is within the +% place field is above 20 % +% --3 Criteria: peaks must be above 50 cm and below 140 cm on the wheel so that +% cells close to reward are not considered as cue tuned + +% +% INPUT +% avgActivity: Array with the length of number of bins on the wheel, +% where each sample is the average activity for all laps. +% A Gaussian smoothing should also have +% been applied to this map (like Mao et al (2017) does). +% +% signalMatrix : Matrix with the no of columns corresponding to no of bins on the wheel +% and no of rows corresponds to no of trials. +% Each sample is the avearge actavity activity for all laps for +% position/ trial bin +% +% OUTPUT +% pfData: A struct containing all related place field information: +% .ifv: mean in-field value +% .field_width: width of place field as number of bins +% .field_peak_ind: index of peak of place field(s) +% .field_start: first index of place field(s) +% .field_end: last index of place field(s) +% .n_peaks: number of peaks in position activity map +% .inFieldIdxs: array containing all indices that are within any field of the position map +% .outOfFieldIdxs: array containing all indices that are considered out of a field +% +% this function was written by Andreas Lande and adapted to cue field +% identification by Ann Christin Garvert + +% Set parameters +% Get initial place field threshold +% Threshold is set as 30% of the difference between the highest and lowest activity of the position activity map. +threshold = min(avgActivity)+(0.3*(max(avgActivity)-min(avgActivity))); +n_bins = 78; +bin_size = sData.behavior.meta.binSize; + +minPeakDist = round(10/bin_size); % Distance needed for two peaks to be considered as individual +minPeakWidth = 1; +maxPeakWidth = 78; +activityRatioCriteria = 3; + +% Check that min peak is not too small + +% Use the MATLAB Signal Processing Toolbox function findpeaks() to localize peaks within the position activity map +[~,locs] = findpeaks(avgActivity,'MinPeakDistance',minPeakDist,'MinPeakHeight',threshold,'WidthReference','halfheight','MinPeakWidth',minPeakWidth,'MaxPeakWidth',maxPeakWidth); + +% Struct for output +pfData = struct(); + +% If no peaks are found, skip, else run analysis +if ~isempty(length(locs) == 0) + + % For each peak detected + for peak = 1:length(locs) + peakIndx = locs(peak); + + % Find all neighbouring values that are above the threshold to obtain + % place field width + pfStart = peakIndx; % start index of place field + pfEnd = peakIndx; % end index of place field + + % Go from peak index backwards to locate beginning of place field + look = 1; + curr_pos = peakIndx; + while look + curr_pos = curr_pos - 1; + if curr_pos < 1 + curr_pos = n_bins; + end + + if avgActivity(curr_pos) > threshold + pfStart = curr_pos; + else + look = 0; + end + end + + % Go forward from peak to find end index of place field + look = 1; + curr_pos = peakIndx; + while look + curr_pos = curr_pos + 1; + if curr_pos > n_bins + curr_pos = 1; + end + + if avgActivity(curr_pos) > threshold + pfEnd = curr_pos; + else + look = 0; + end + end + + % Get all in-field and out-of-field indices + inFieldIdxs = []; + outOfFieldIdxs = []; + if pfEnd < pfStart + inFieldIdxs = [inFieldIdxs, 1:pfEnd]; + inFieldIdxs = [inFieldIdxs, pfStart:length(avgActivity)]; + outOfFieldIdxs = [pfEnd+1:pfStart-1]; + else + inFieldIdxs = [inFieldIdxs, pfStart:pfEnd]; + outOfFieldIdxs = [1:pfStart-1,pfEnd+1:length(avgActivity)]; + end + + inFieldValue = nanmean(avgActivity(inFieldIdxs)); +% out_of_field_value = nanmean(pos_activity_map(outOfFieldIdxs)); + + % Save place field information + pfData.inFieldValue(peak) = inFieldValue; + pfData.fieldWidth(peak) = length(inFieldIdxs); + pfData.fieldPeakIdx(peak) = locs(peak); + pfData.fieldStart{peak} = pfStart; + pfData.fieldEnd{peak} = pfEnd; + pfData.nPeaks = length(locs); + pfData.inFieldIdxs{peak} = inFieldIdxs; + pfData.outOfFieldIdxs{peak} = outOfFieldIdxs; + + + end + + % For each field, obtain out-of-field/in-field ratio + totalOutOfFieldIdx = 1:n_bins; % Start with the assumption that all bins are out of field + totalInFieldIdxs = [pfData.inFieldIdxs{:}]; + pfData.totalInFieldIdxs = totalInFieldIdxs; + totalOutOfFieldIdx = totalOutOfFieldIdx(~ismember(totalOutOfFieldIdx,totalInFieldIdxs)); + pfData.totalOutOfFieldIdx = totalOutOfFieldIdx; + pfData.outOfFieldValue(peak) = mean(avgActivity(pfData.totalOutOfFieldIdx)); + pfData.inFieldValue(peak) = mean(avgActivity(pfData.totalInFieldIdxs)); + pfData.activityRatioInVsOut(peak) = pfData.inFieldValue(peak)/pfData.outOfFieldValue(peak); + + + % calculate for each trials in which activity peak is within the place field + pfData.isAboveThresholdTrial(peak) = 0; + for lap = 1:numel(sData.trials.trialStart)-1 + pfData.inFieldValue_trial(peak,lap) = mean(signalMatrix(lap,pfData.inFieldIdxs{peak})); + pfData.outOfFieldValue_trial(peak,lap) = mean(signalMatrix(lap,pfData.totalOutOfFieldIdx)); + + if pfData.inFieldValue_trial(peak,lap)/pfData.outOfFieldValue_trial(peak,lap) > activityRatioCriteria + pfData.isAboveThresholdTrial(peak) = pfData.isAboveThresholdTrial(peak) +1; + end + end + + % correct for fields with multiple peaks + [startVals, startIdx, startValsOcc]= unique([pfData.fieldStart{:}]); + deleteIdx = []; + if length(startVals) < length([pfData.fieldStart{:}]) + for i = 1:length(startIdx) + noPeak = find(startValsOcc == startIdx(i)); + if length(noPeak) > 1 + deleteIdx = [deleteIdx; noPeak(2:end)]; + + [~,maxPeakIdx] = max(avgActivity(noPeak)); + + keepIdx = pfData.fieldPeakIdx(noPeak(maxPeakIdx)); + pfData.fieldPeakIdx(noPeak(1)) = keepIdx; + + + pfData.nPeaks = pfData.nPeaks-(length(noPeak)-1); + end + end + + pfData.inFieldValue(deleteIdx) = []; + pfData.fieldWidth(deleteIdx) = []; + pfData.fieldStart(deleteIdx) = []; + pfData.fieldEnd(deleteIdx) = []; + pfData.inFieldIdxs(deleteIdx) = []; + pfData.outOfFieldIdxs(deleteIdx) = []; + pfData.fieldPeakIdx(deleteIdx) = []; + + end + + + % Based on the activity ratio measure obtained above peaks must fulfill the + % following: + % --1 Criteria: the activityRatio of in field activity vs out of field activity is above 3 (ActivityRatioCriteria = 3) % This is the ratio used by Mao et al (2017) + % --2 Criteria: The percentage of laps where the activity peak is within the place field is above 30 % + % --3 Criteria: peaks must be above 50 cm and below 140 cm on the wheel so that cells close to reward are not considered as cue tuned + + for peak = 1:pfData.nPeaks + pfData.isActive(peak) = 1; + % --1 Criteria: + if pfData.activityRatioInVsOut(peak) < activityRatioCriteria + pfData.isActive(peak) = 0; end + + % --2 Criteria: + if pfData.isAboveThresholdTrial(peak)/(numel(sData.trials.trialStart)-1) < 0.2 + pfData.isActive(peak) = 0; + end + + % --3 Criteria + if pfData.fieldPeakIdx(peak) < 25 || pfData.fieldPeakIdx(peak) > 70 + pfData.isActive(peak) = 0; + end + + + if pfData.isActive(peak) == 0 + + pfData.inFieldValue(peak)= []; + pfData.fieldWidth(peak) = []; + pfData.fieldStart(peak) = []; + pfData.fieldEnd(peak) = []; + pfData.inFieldIdxs(peak) = []; + pfData.outOfFieldIdxs(peak) = []; + pfData.fieldPeakIdx(peak) = []; + pfData.isAboveThresholdTrial(peak)= []; + + pfData.nPeaks = pfData.nPeaks-1; + pfData.isActive(peak) = []; + end + + end + +end \ No newline at end of file diff --git a/analysis/cellClassification/fit_model.m b/analysis/cellClassification/fit_model.m new file mode 100644 index 0000000..6d56fd6 --- /dev/null +++ b/analysis/cellClassification/fit_model.m @@ -0,0 +1,60 @@ +function [testFit,trainFit,param_mean] = fitModel(stateMatrix,signal,beta) + +% make signal and stateMatrix same size: +stateMatrix = stateMatrix(1:min([size(stateMatrix,1),size(signal,2)]),:); +signal = signal(1:min([size(stateMatrix,1),size(signal,2)])); + +%% Initialize matrices and section the data for k-fold cross-validation +numFolds = 5; +[~,numCol] = size(stateMatrix); +sections = numFolds*5; + +% divide the data up into 5*num_folds pieces +edges = round(linspace(1,numel(signal)+1,sections+1)); + +% initialize matrices +testFit = nan(numFolds,3); % llh increase, log_llh_test, log_llh_rand_test +trainFit = nan(numFolds,3); % llh increase, log_llh_train, log_llh_rand_train +paramMat = nan(numFolds,numCol); + +%% perform 5-fold cross validation +for k = 1:numFolds + + %% prepare test and train data sets + testInd = [edges(k):edges(k+1)-1 edges(k+numFolds):edges(k+numFolds+1)-1 ... + edges(k+2*numFolds):edges(k+2*numFolds+1)-1 edges(k+3*numFolds):edges(k+3*numFolds+1)-1 ... + edges(k+4*numFolds):edges(k+4*numFolds+1)-1] ; + + % test data + testSignal = signal(testInd)'; + testMatrix = stateMatrix(testInd,:); + + % training data + trainInd = setdiff(1:numel(signal),testInd); + trainSignal = signal(trainInd)'; + trainMatrix = stateMatrix(trainInd,:); + + %% fit the parameters + opts = optimset('Gradobj','off','Hessian','off','Display','off'); + + if k == 1 + init_param = 1e-3*randn(numCol, 1); + else + init_param = param; + end + + init_param(isinf(init_param)) = 0; init_param(isnan(init_param)) = 0; + param = fminunc(@(param) lnPoissonModel(param,beta,trainMatrix,trainSignal),init_param,opts); + + + % calculate log-likelihood increase from "mean firing rate model" + testFit(k,:) = calculateLLHIncrease(testSignal,testMatrix,param); + trainFit(k,:) = calculateLLHIncrease(trainSignal,trainMatrix,param); + paramMat(k,:) = param; + +end + +param_mean = nanmean(paramMat); + +return +end \ No newline at end of file diff --git a/analysis/cellClassification/main_cell_classification.m b/analysis/cellClassification/main_cell_classification.m new file mode 100644 index 0000000..9c2eeb8 --- /dev/null +++ b/analysis/cellClassification/main_cell_classification.m @@ -0,0 +1,59 @@ + + +%% Input structure (required) +% signal : array of deconvolved signal for neurons x length of number time bins +% covariate struct with fields: +% {'type' } % name of covariate +% {'predictor' } % variable of this covariate for each time bin +% {'no bins' } % categorical bins that will be created +% {'beta' } % this is the smoothing parameter +% {'circular' } % is the max bin equal to min? +% + +%% Output structure +% results with fiels: +% {'testFit' } : log likelihood for test data: folds x covariate +% {'trainFit' } : log likelihood for train data: folds x covariate +% {'params' } : fitted model parameters : covariate x no bins +% {'selected_model' } : integer corresponding to cov struct + + +% llh increase, log_llh_test, log_llh_rand_test + +% other variables +noCov = length(cov); % number of covariates tested +numFold = 5; % attention this is also hard coded in the function fit_model + +% prepare Signal +numNeuron = size(signal,1); + +% create state matrix : +for j = 1:noCov + covGrid{j} = stateMatrix(cov(j)); % position +end + +% calculate LLH +results = struct(); + +for i = 1:numNeuron + + for j = 1:noCov + [results(i).testFit{j},results(i).trainFit{j},results(i).param{j}] = fit_model(covGrid{j},signal(i,:), cov(j).beta); + end + + results(i) = compareModelPerformance(results(i),numFold,noCov); + [isCue, results(i).cueData] = checkForCueTuning(sData,signal(i,:)); + + % if is cue cell, the + if isCue == 1 + results(i).selected_model = noCov+1; + end + + +end + + + + + + \ No newline at end of file diff --git a/analysis/cellClassification/run_example.m b/analysis/cellClassification/run_example.m new file mode 100644 index 0000000..0e768ee --- /dev/null +++ b/analysis/cellClassification/run_example.m @@ -0,0 +1,55 @@ + +%% Clear the workspace and load the data + +% clear all; close all; clc +% load sData with data + +%% Input structure (required) +% signal : No. neurons x No. Time points -matrix +% make sure that the signal and all predictors have the same length! + +signal = double(sData.imdata.roiSignals(2).dff); + +% position +cov(1).type = 'position'; +cov(1).predictor = sData.behavior.wheelPosDs; +cov(1).bins = 40; +cov(1).beta = 0.01; +cov(1).limits = [0 157]; +cov(1).circular = 'true'; + +% speed +cov(2).type = 'speed'; +cov(2).predictor = sData.behavior.runSpeedDs; +cov(2).bins = 30; +cov(2).beta = 0.01; +cov(2).limits = [0 60]; +cov(2).circular = 'false'; + +% acceleration +cov(1).type = 'acceleration'; +dt = 1/31; +cov(3).predictor = diff(sData.behavior.runSpeedDs)./dt; +cov(3).bins = 40; +cov(3).beta = 0; +cov(3).limits = [-200 200]; +cov(3).circular = 'false'; + +% lick rate +cov(4).type = 'lick_rate'; +cov(4).lickLengthThreshold = 15; +cov(4).lickFreqThreshold = 15; +cov(4).predictor = extractLickFrequency(sData, 10000,cov(4).lickLengthThreshold, cov(4).lickFreqThreshold); +cov(4).bins = 20; +cov(4).beta = 0; +cov(4).limits = [0 10]; +cov(4).circular = 'false'; + + +main_cell_classification + + + + + + \ No newline at end of file diff --git a/analysis/cellClassification/splitInTrials.m b/analysis/cellClassification/splitInTrials.m new file mode 100644 index 0000000..aa37174 --- /dev/null +++ b/analysis/cellClassification/splitInTrials.m @@ -0,0 +1,15 @@ +function signalMatrix = splitInTrials(signal, sData) +% this script organises into a matrix in which each row corresponds to a +% trial and each column to a position on the track (2cm binning) + + signalMatrix = zeros(numel(sData.trials.trialStart)-1,floor(sData.stats.wheelCircum/2)); + for i = 1:numel(sData.trials.trialStart)-1 + isInTrial = sData.trials.trialLength == i; + trialSignal = signal(isInTrial); + + for jj = 1:floor(sData.stats.wheelCircum/2) + idxInTrialBin = sData.behavior.wheelPosDsBinned(isInTrial) == jj; + signalMatrix(i,jj) = nanmean(trialSignal(idxInTrialBin)); + end + end +end \ No newline at end of file diff --git a/analysis/cellClassification/stateMatrix.m b/analysis/cellClassification/stateMatrix.m new file mode 100644 index 0000000..8b84656 --- /dev/null +++ b/analysis/cellClassification/stateMatrix.m @@ -0,0 +1,37 @@ +function covGrid = stateMatrix(covariat)%,limits,binNo,iscircular) +% changes covariate into categorical numbers and creates state matrix +% input: +% covariate: eg: position, speed, acceleration or lick rate +% binNo # states for covariate +% iscircular: boolean that declares wether covariate is circular (as +% position i +% cov(1).type = 'position'; +% cov(1).predictor = sData.behavior.wheelPosDs; +% cov(1).bins = 40; +% cov(1).beta = 0.01; +% cov(1).limits = [0 157]; +% cov(1).circular = 'true'; + +% create state matrix +binSize = (covariat.limits(2)-covariat.limits(1))/covariat.bins; + +covariatBinned = zeros(length(covariat.predictor),1); + +if covariat.circular + binGaps = binSize/2:binSize:covariat.limits(2)-binSize/2; +else + binGaps = covariat.limits(1):binSize:covariat.limits(2); +end + +for t = 1:length(covariat.predictor) + [~,ind] = min(abs(binGaps-covariat.predictor(t))); + covariatBinned (t) = ind; +end + +covGrid = zeros(length(covariatBinned),length(binGaps)+1); +for l = 1: length(covariatBinned) + covGrid(l,covariatBinned(l)+1) = 1; +end +covGrid(:,1) = 1; + +end \ No newline at end of file diff --git a/analysis/decoding/calcMaxLikelihoodPos.m b/analysis/decoding/calcMaxLikelihoodPos.m new file mode 100644 index 0000000..f89f364 --- /dev/null +++ b/analysis/decoding/calcMaxLikelihoodPos.m @@ -0,0 +1,16 @@ + +function predPos = calcMaxLikelihoodPos(prior,likelihood, binGaps) + + predPos = nan(size(likelihood,1),1); + posterior_NotNorm = nan(size(likelihood)); + posterior = nan(size(likelihood)); + + for t = 1:size(likelihood_train,1) + posterior_NotNorm(t,:) = likelihood(t,:).*prior; + posterior(t,:) = posterior_NotNorm(t,:)./nansum(posterior_NotNorm(t,:)); + + [~, maxInd] = nanmax(posterior(t,:)); + predPos(t) = binGaps(maxInd); + end + +end \ No newline at end of file diff --git a/analysis/decoding/calculateDecodingError.m b/analysis/decoding/calculateDecodingError.m new file mode 100644 index 0000000..55079af --- /dev/null +++ b/analysis/decoding/calculateDecodingError.m @@ -0,0 +1,16 @@ +function rmse = calculateDecodingError(predPos, realPos) + + maxval = max(realPos); + predError = zeros(length(realPos),1); + for i = 1:length(realPos) + predError(i) = abs(realPos(i)-predPos(i)); + if predError(i) > maxval/2 && realPos(i) > predPos(i) + predError(i) = abs(maxval-realPos(i)+predPos(i)); + elseif predError(i) > maxval/2 && realPos(i) < predPos(i) + predError(i) = abs(maxval-predPos(i)+realPos(i)); + end + end + + rmse = mean(predError); + +end \ No newline at end of file diff --git a/analysis/decoding/decode_main.m b/analysis/decoding/decode_main.m new file mode 100644 index 0000000..747e980 --- /dev/null +++ b/analysis/decoding/decode_main.m @@ -0,0 +1,114 @@ +close all + +file = 'm5094-20200522-02'; +filePath = ''; +sData = load(fullfile(filePath,file,'/final_sData.mat')); + +savePath = ''; +if ~exist(savePath, 'dir'); mkdir(savePath); end + + +%% set variables +numFolds = 5; % number of crossvalidation folds + +% hyperparameters that were optimized beforehand: +tempBinSize = 3; % temporal bin sized +spatialbinSize = 5.233; % so that we have 30 evenly sized spatial bins +beta = 1; % roughness penalty +speedThreshold = 1; + +% initialise signal and position +deconvSignalTemp = sData.imdata.roiSignals(2).deconv; +pos = sData.behavior.wheelPosDs; + +% create time bins +timeBins = 1:tempBinSize:length(pos); +for t = 1:length(timeBins)-1 + position(t) = nanmean(pos(timeBins(t):timeBins(t+1)-1)); + signal(:,t) = nanmean(deconvSignalTemp(:,timeBins(t):timeBins(t+1)-1),2); + speed(t) = nanmean(sData.behavior.runSpeedDs(timeBins(t):timeBins(t+1)-1)); +end + + +% make sure all parameters have the same length +minLength = min([length(position), size(signal,2), length(speed), length(sData.trials.trialLength)]); +position = position(1:minLength); +signal = signal(:,1:minLength); +speed = speed(1:minLength); + + +binGaps = spatialbinSize/2:spatialbinSize:max(position); +bins = 1:length(binGaps); +posBinned = zeros(length(position),1); + +for t = 1:length(position) + [~,ind] = min(abs(binGaps-position(t))); + posBinned(t) = ind; +end + + +% delete indices where animal is too slow +signal(:,(runSpeed < speedThreshold)) = []; +posBinned(runSpeed < speedThreshold) = []; + +% divide the data up into 5*num_folds pieces + +sections = numFolds*5; +edges = round(linspace(1,size(signal,2)+1,sections+1)); + +d_data = []; + +for k = 1:numFolds + testIdx = [edges(k):edges(k+1)-1 edges(k+numFolds):edges(k+numFolds+1)-1 ... + edges(k+2*numFolds):edges(k+2*numFolds+1)-1 edges(k+3*numFolds):edges(k+3*numFolds+1)-1 ... + edges(k+4*numFolds):edges(k+4*numFolds+1)-1]; + + trainIdx = setdiff(1:size(signal,2),testIdx); + + d_data.iter{k}.trainSignal = signal(:,trainIdx); + d_data.iter{k}.testSignal = signal(:, testIdx); + + d_data.iter{k}.realPos_train = posBinned(trainIdx); + d_data.iter{k}.realPos_test = posBinned(testIdx); + + %% find cells without signal and delete + d_data.iter{k}.cellIdx = 1:size(signal,1); + d_data.iter{k} = findEmptyROIs(d_data.iter{k}); + + %% calculate encoding models for + d_data.iter{k}.lambda = fitModelDec(d_data.iter{k}.trainSignal, posBinned(trainIdx),beta); + + % calculate likelihood: P(signal|position) + likelihood_test = likelihood(d_data.iter{k}.testSignal, d_data.iter{k}.lambda); + likelihood_train = likelihood(d_data.iter{k}.trainSignal, d_data.iter{k}.lambda); + + + %% calculate posterior: P(position| signal) = P(signal|position)*P(position) + % calculate prior P(position) --> flat prior, assume every position is equally likely + % predicted position is the position with maximum likelihood + prior = ones(1,size(likelihood_train,2)).*1/size(likelihood_train,2); + + d_data.iter{k}.predPos_train = calcMaxLikelihoodPos(prior,likelihood_train, binGaps); + d_data.iter{k}.predPos_test = calcMaxLikelihoodPos(prior,likelihood_test, binGaps); + + d_data.rmse_test(k) = calculateDecodingError(d_data.iter{k}.predPos_test, d_data.iter{k}.realPos_test); + d_data.rmse_train(k) = calculateDecodingError(d_data.iter{k}.predPos_train, d_data.iter{k}.realPos_train); + + +end + + % choose fold with median performance in test data: +d_data.medianRMSE_test = nanmedian(d_data.rmse_test(:)); +[~, k] = min(abs(d_data.rmse_test(:)-median(d_data.rmse_test(:)))); + +% save corresponding train rmse's +d_data.medianRMSE_train = d_data.rmse_train(k); + +d_data.meanRMSE_test = nanmean(rmse_test); +d_data.meanRMSE_train = nanmean(rmse_train); + +d_data.medianTestFold = k; +d_data.lambda = d_data.iter{k}.lambda; + + +save(fullfile(savePath, '/d_data.mat'), 'd_data','-v7.3'); diff --git a/analysis/decoding/findEmptyROIs.m b/analysis/decoding/findEmptyROIs.m new file mode 100644 index 0000000..86466fa --- /dev/null +++ b/analysis/decoding/findEmptyROIs.m @@ -0,0 +1,15 @@ +function data = findEmptyROIs(data) + + deleteCell = []; + for n = 1:length(data.cellIdx) + if isempty(find(data.trainSignal(n,:) ~= 0, 1))|| isempty(find(data.testSignal(n,:)~= 0, 1))||... + length(find(isnan(data.trainSignal(n,:)) == 1)) == length(find(isnan(data.trainSignal(n,:)) ... + || length(find(isnan(data.testSignal(n,:)) == 1)) == length(find(isnan(data.testSignal(n,:)) + deleteCell(end+1) = n; + end + end + + data.deconvTrain(deleteCell,:) = []; data.deconvTest(deleteCell,:) = []; + data.cellIdx(deleteCell) = []; + +end \ No newline at end of file diff --git a/analysis/decoding/fitModelDec.m b/analysis/decoding/fitModelDec.m new file mode 100644 index 0000000..7d7ad6e --- /dev/null +++ b/analysis/decoding/fitModelDec.m @@ -0,0 +1,24 @@ +function lambda = fitModelDec(signal, pos,beta) + +% create state matrix +X = zeros(length(pos),length(unique(pos)+1)); +for j = 1:length(pos); X(j,pos(j)+1) = 1; end; X(:,1) = 1; + +% initialize firing rate lambda +lambda = zeros(size(signal,1),max(pos)); +% fit parameters +for i = 1:size(signal,1) + + opts = optimset('Gradobj','off','Hessian','off','Display','off'); + + init_param = (X'*signal(i,:)')/sum(signal(i,:)); + + init_param(isinf(init_param)) = 0; init_param(isnan(init_param)) = 0; + + lossfun = @(params) lnPoissonModel(params,beta, X, signal(i,:)'); + param = fminunc(lossfun,init_param,opts); + + lambda = exp(X * param'); +end + +end \ No newline at end of file diff --git a/analysis/decoding/likelihood.m b/analysis/decoding/likelihood.m new file mode 100644 index 0000000..8bc5a45 --- /dev/null +++ b/analysis/decoding/likelihood.m @@ -0,0 +1,28 @@ + +function likelihood = likelihood(signal,lambda) +% lambda1(n,:) = b(n)+g(n)*exp(k(n)*cos(x-mu(n))); +% calculate P(spike(i)| time(i)) + +prob = zeros(size(lambda)); +likelihood = zeros(size(signal,2),size(lambda,2)); + +for t = 1: size(deconv,2) + % select signal for one time bin + frTemp = repmat(signal(:,t),1,size(lambda,2)); + + for j = 1:size(frTemp,2) + for i = 1:size(frTemp,1) + %prob. of the time points signals (separate per neuron) given position + prob(i,j) = exp(frTemp(i,j) * log(lambda(i,j)) - lambda(i,j) - gammaln(frTemp(i,j)+1)); + end + end + + % probability signal(each time point) given position, multiply over + % neurons -> here we assume statistical independence + likelihood(t,:) = prod(prob); %Prob. of all signals (all neurons) given pos +end + +% normalisation +likelihood = likelihood./nansum(likelihood,2); + +end \ No newline at end of file diff --git a/analysis/glm/lnPoissonModel.m b/analysis/glm/lnPoissonModel.m new file mode 100644 index 0000000..bddf557 --- /dev/null +++ b/analysis/glm/lnPoissonModel.m @@ -0,0 +1,16 @@ +function logL = lnPoissonModel(param,b,X,Y) +% X: covariate state matrix +% Y: signal + +nonlinear = @(x) exp(x) ; +lambda = nonlinear(X * param); % firing rate + +% compute smoothing of parameters +smoothVal = smoothParams(param,b); + +% compute f, the gradient, and the hessian +logL = sum(lambda-Y(:).*log(lambda)) + smoothVal; + +end + + \ No newline at end of file diff --git a/analysis/glm/smoothParams.m b/analysis/glm/smoothParams.m new file mode 100644 index 0000000..4cb82d7 --- /dev/null +++ b/analysis/glm/smoothParams.m @@ -0,0 +1,7 @@ +function smoothMatrix = smoothParams(param,beta) + % smoothing penalty + numParam = numel(param); + D1 = spdiags(ones(numParam,1)*[-1 1],0:1,numParam-1,numParam); + DD1 = D1'*D1; + smoothMatrix = beta*0.5*param'*DD1*param; +end \ No newline at end of file diff --git a/analysis/skaggs/calculate_skaggs.m b/analysis/skaggs/calculate_skaggs.m new file mode 100644 index 0000000..901c145 --- /dev/null +++ b/analysis/skaggs/calculate_skaggs.m @@ -0,0 +1,113 @@ +function [SI, sig, pVal ,specificity] = calculate_skaggs(occupancy_map_2D, activity_map_2D) +%% CALC_SPATIALINFORMATION +% Calculates spatial information according to the method used by Skaggs et al (1993) +% +% INPUT +% occupancy_map: A matrix where each row corresponds to a lap of +% running and each column is a binned position on the wheel. Each +% element contains the number of samples the animal was in the bin. +% activity_map_2D: A matrix where each row corresponds to a lap of +% running and each column is a binned position on the wheel. Each +% element contains summed activity of a signal (fex deconvolved) +% when the animal was in the bin. +% +% OUTPUT +% SI: The spatial information as bits. +% sig: Outputs result of a significance test. 1 if significant, 0 else. +% For the SI to be significant, it has to be higher than the 95th +% percentile of SI calculated for 1000 shuffled datasets. + + +%% Parameters +sig= []; pVal= []; +% bin_size = 1; %cm +n_bins = size(activity_map_2D,2); % The track number of bins + +% Down project maps from 2D to 1D +occupancy_map_1D = sum(occupancy_map_2D,1); +% activity_map_1D = sum(activity_map_2D,1); + +%% Calculate spatial information +% Create a position corrected activity map +% position_activity_map_1D = activity_map_1D./occupancy_map_1D; +mean_firing_of_bin_1D = nanmean(activity_map_2D); + +% Smooth position activity map using a Gaussian window if deconv is used +% position_activity_map_1D = smoothdata(position_activity_map_1D,'gaussian',7); % position_activity_map_1D; % + +% Create occupancy proability array +occupancy_probability = zeros(size(occupancy_map_1D)); +for p = 1:length(occupancy_map_1D) + occupancy_probability(p) = occupancy_map_1D(p)/nansum(occupancy_map_1D); +end + +SI = NaN(n_bins,1); +for bin = 1:n_bins + SI(bin) = occupancy_probability(bin) * (mean_firing_of_bin_1D(bin)/nanmean(mean_firing_of_bin_1D)) * ... + log2(mean_firing_of_bin_1D(bin)/mean(mean_firing_of_bin_1D)); +end + + +% Sum the spatial information across all bins to get the spatial +% information in bits +if nansum(SI) < 0 + error('Summed SI is below zero'); +end + +SI = nansum(SI)./0.14; +% information rate divided by overall mean firing rate leads to information +% per spike/specificity +specificity = SI/mean(mean_firing_of_bin_1D); + +%% Run significance test by shuffling data +% n_shuffles = 1000; +% n_laps = size(occupancy_map_2D,1); +% +% % Initiate shuffled matrix +% shuffled_SI = []; +% +% % Make n_shuffles different shuffles +% for i = 1:n_shuffles +% colShifts = randi(n_bins, 1, n_laps); +% activity_shuffled = zeros(n_laps,n_bins); +% occup_shuffled = zeros(n_laps,n_bins); +% +% for l = 1:n_laps +% activity_shuffled(l,:) = circshift(activity_map_2D(l,:),colShifts(l)); +% occup_shuffled(l,:) = circshift(occupancy_map_2D(l,:), colShifts(l)); +% end +% +% % Use shuffled activity to estimate new SI +% % Downproject into 1D +% activity_shuffled = sum(activity_shuffled,1); +% +% occup_shuffled_map_1D = sum(occup_shuffled,1); +% +% % Create a position corrected activity map +% position_activity_map_1D = activity_shuffled./occup_shuffled_map_1D; +% +% % Smooth position activity map using a Gaussian window with std 4.5 cm +% % given a bin size of 1.5 cm +% %position_activity_map_1D = smoothdata(position_activity_map_1D,'gaussian',3); +% +% SI_shuff = []; +% for bin = 1:n_bins +% SI_shuff(bin) = occupancy_probability(bin) * (position_activity_map_1D(bin)/mean(position_activity_map_1D)) * ... +% log2(position_activity_map_1D(bin)/mean(position_activity_map_1D)); +% end +% +% % Sum the spatial information across all bins to get the spatial +% % information in bits +% shuffled_SI(i) = nansum(SI_shuff); +% +% end +% +% pVal = length(find(shuffled_SI > SI))/1000; +% % Determine significance by using 95th percentile +% if SI > prctile(shuffled_SI,99) +% sig = 1; +% else +% sig = 0; +% end + +end \ No newline at end of file