diff --git a/Code/ephysMain.m b/Code/ephysMain.m index 932399d..d6cff1c 100644 --- a/Code/ephysMain.m +++ b/Code/ephysMain.m @@ -1451,19 +1451,26 @@ xtickangle(45); %% region comparisons of cue, lick, and both cells using generalized linear mixed effects model figure; +pcutoff=0.05/36; regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; -reggrps=[1:3]; +reggrps=[2:4]; +neuronGroupAdj = neuronGroup + 1; +neuronGroupAdj(neuronGroupAdj==5)=1; allReg=unique(neuronRegionOlf); regInd=NaN(length(regions),1); +order=[9 10 5 2 3 1 7 12 8 11 6 13 4]; +allReg=allReg(order); for reg=1:length(regions) regInd(reg,1)=find(ismember(allReg,regions(reg))); end + %fit all together, allowing random effect of session to have more power catcolors={[0 0.4 0.9],[0.9 0.2 0],[0.6 0 0.6]}; titles={'cues','licks','both'}; cats=[1 3 2]; +pVals=NaN(36,length(cats)); for ct=1:3 %region barData=zeros(2,5); @@ -1477,9 +1484,40 @@ barData(1,reg)=sum(sel®sel)/sum(regsel); end tbl = table(categorical(neuronSession(regionnumber>0)), categorical(regionnumber(regionnumber>0)), sel(regionnumber>0), 'VariableNames', {'session', 'region','val'}); - lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial'); + lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial','dummyvarcoding','reference'); [ypred,ypredCI]=predict(lm,'conditional',false); + Hall=[]; +% Htemp = [1 0 0 0 0 0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0 0 0 0 0 0]; +% for reg1 = 1:length(regions) +% H=Htemp; +% H(2,reg1+4)=1; +% Hall=cat(1,Hall,H); +% end + Htemp = [0 0 0 0 0 0 0 0 0 0 0 0 0]; + for reg1 = 1:length(regions)-1 + for reg2 = reg1+1:length(regions) + H=Htemp; + H(reg1+4)=1; + H(reg2+4)=-1; + Hall=cat(1,Hall,H); + end + end + greater=zeros(length(regions)); + for reg1 = 1:length(regions) + for reg2 = 1:length(regions) + H=Htemp; + H(reg1+4)=1; + H(reg2+4)=-1; + greater(reg1,reg2)=coefTest(lm,H); + end + end + + %Hall=unique(Hall,'rows'); + + for p=1:length(pVals) + pVals(p,ct) = coefTest(lm,Hall(p,:)); + end barData=barData(:,regInd); b=bar(barData','stacked'); %b(3).FaceColor=[0.4 0 0.5]; %odor @@ -1517,16 +1555,17 @@ customMap=cat(1,[1 1 1],customMap); %by region - greater=zeros(length(regions)); + %greater=zeros(length(regions)); frac=zeros(length(regions)); for r1=1:length(regions) for r2=1:length(regions) frac(r1,r2)=barData(1,r1)-barData(1,r2); - greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); + %greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end s=subplot(4,3,3+ct); - imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); + + imagesc(frac .* (greater0), [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); cb=colorbar; @@ -1549,7 +1588,7 @@ sel=ismember(category,cats(ct)); regionnumber=zeros(totalNeurons,1); for reg=1:4 - regsel=ismember(neuronGroup,reg); + regsel=ismember(neuronGroupAdj,reg); regionnumber(regsel)=reg; barData(1,reg)=sum(sel®sel)/sum(regsel); end @@ -1569,7 +1608,7 @@ lowerCI=[]; estmean=[]; for reg=1:length(reggrps) - regsel=find(ismember(neuronGroup,reggrps(reg)),1); + regsel=find(ismember(neuronGroupAdj,reggrps(reg)),1); estmean(1,reg)=ypred(regsel); upperCI(1,reg)=ypredCI(regsel,2); lowerCI(1,reg)=ypredCI(regsel,1); @@ -1580,11 +1619,21 @@ xlim([0.5 7+0.5]); %legend('Cues','Cues,Licks','Licks','Licks,Reward','Reward','Cues,Reward','Cues,Licks,Reward'); xticks(1:length(reggrps)); - xticklabels(regionGroupNames(reggrps)); + xticklabels(regionGroupNames(reggrps-1)); xtickangle(45); ylim([0 0.5]); - + Htemp = [0 0 0 0]; + greater=zeros(length(reggrps)); + for reg1 = 1:length(reggrps) + for reg2 = 1:length(reggrps) + H=Htemp; + H(reg1+1)=1; + H(reg2+1)=-1; + greater(reg1,reg2)=coefTest(lm,H); + end + end + %make heatmap target=catcolors{ct}; customMap=NaN(100,3); @@ -1594,25 +1643,25 @@ customMap=cat(1,[1 1 1],customMap); %by region - greater=zeros(3); + %greater=zeros(3); frac=zeros(3); for r1=1:length(reggrps) for r2=1:length(reggrps) frac(r1,r2)=barData(1,r1)-barData(1,r2); - greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); + %greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end s=subplot(4,3,9+ct); - imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); + imagesc(frac .* (greater<0.05/3 & frac>0), [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); cb=colorbar; set(cb,'ticks',[0:0.1:0.3]); xticks(1:length(reggrps)); - xticklabels(regionGroupNames(reggrps)); + xticklabels(regionGroupNames(reggrps-1)); xtickangle(45); yticks(1:length(reggrps)); - yticklabels(regionGroupNames(reggrps)); + yticklabels(regionGroupNames(reggrps-1)); title(titles{ct}); if ct==1 ylabel('increase in this region...'); @@ -1621,6 +1670,8 @@ end + + %% visualize GLM fitting %1380 @@ -2224,7 +2275,207 @@ end end +%% trying different GLM cutoffs +pcutoff=0.05/36; +figure; + +map=redblue(256); + +regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; +neuronRegionOlf=neuronRegionAdj; +neuronRegionOlf(ismember(neuronRegionAdj,{'EPd','PIR'}))={'OLF'}; +allReg=unique(neuronRegionOlf); +order=[9 10 5 2 3 1 7 12 8 11 6 13 4]; +allReg=allReg(order); +for reg=1:length(regions) + regInd(reg,1)=find(ismember(allReg,regions(reg))); +end + + +glmcutoffs = [0.005 0.01;0.01 0.02;0.02 0.03;0.03 0.05;0.05 1]; +sels={}; +for g=1:length(glmcutoffs) + cueSel = improvement(:,2)>glmcutoffs(g,1) & varExp(:,1)>glmcutoffs(g,1); + + subplot(6,5,g); + cueSel=cueSel & improvement(:,3)<=glmcutoffs(g,1) & improvement(:,4)<=glmcutoffs(g,1); + hold on + barData=zeros(2,length(regions)); + for reg=1:length(regions) + regsel=ismember(neuronRegionAdj,regions(reg)); + barData(1,reg)=sum(cueSel®sel)/sum(regsel); + end + + b=bar(barData','stacked'); + b(1).FaceColor=[0 0.4 0.9]; %cue + b(2).FaceColor=[0.8 0.8 0.9]; + + xlim([0.5 length(regions)+0.5]); + ylim([0 0.55]); + xticks(1:length(regions)); + xticklabels(regions); + xtickangle(45); + + title(sprintf('>%g%% cue variance (%g%%)',100*(glmcutoffs(g)),round(100*sum(cueSel)/length(cueSel),2,'significant'))); + if g==1 ylabel('fraction of neurons'); end + sel=cueSel; + +if sum(sel)>50 + + regionnumber=zeros(totalNeurons,1); + for reg=1:length(allReg) + regsel=ismember(neuronRegionOlf,allReg(reg)); + regionnumber(regsel)=reg; + end + tbl = table(categorical(neuronSession(regionnumber>0)), categorical(regionnumber(regionnumber>0)), sel(regionnumber>0), 'VariableNames', {'session', 'region','val'}); + lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial','dummyvarcoding','reference'); + [ypred,ypredCI]=predict(lm,'conditional',false); + + greater=zeros(length(regions)); + for reg1 = 1:length(regions) + for reg2 = 1:length(regions) + H=Htemp; + H(reg1+4)=1; + H(reg2+4)=-1; + greater(reg1,reg2)=coefTest(lm,H); + end + end + + upperCI=[]; + lowerCI=[]; + estmean=[]; + for reg=1:length(regions) + regsel=find(ismember(neuronRegionOlf,regions(reg)),1); + estmean(1,reg)=ypred(regsel); + upperCI(1,reg)=ypredCI(regsel,2); + lowerCI(1,reg)=ypredCI(regsel,1); + end + + errorbar(1:length(regions),estmean(1,:),estmean(1,:)-lowerCI(1,:),upperCI(1,:)-estmean(1,:),'o','linewidth',0.75,'color','k'); + + %make heatmap + target=[0 0.4 0.9]; + customMap=NaN(100,3); + for i=1:length(customMap) + customMap(i,:)=target*i/(length(customMap)); + end + customMap=cat(1,[1 1 1],customMap); + + %by region + %greater=zeros(length(regions)); + frac=zeros(length(regions)); + for r1=1:length(regions) + for r2=1:length(regions) + frac(r1,r2)=barData(1,r1)-barData(1,r2); + %greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); + end + end + s=subplot(6,5,5+g); + imagesc(frac .* (greater0), [0 ceil(max(frac,[],'all')*20)/20]); + + colormap(s,customMap); + cb=colorbar; + set(cb,'ticks',[0:0.1:0.3]); + xticks(1:length(regions)); + xticklabels(regions); + xtickangle(45); + yticks(1:length(regions)); + yticklabels(regions); + if g==1 + ylabel('increase in this region...'); + xlabel('over this region'); + end + + %plot heatmaps of cue neurons in each group + +end + +sels{g} = sel; +end + + +%heatmaps +cspActivity=dffPSTH3{1,1}; %get the firing rates of neurons of interest +cueWindow=[0 1.5]; +cueBins=binTimes>=cueWindow(1) & binTimes<=cueWindow(2); + + +colors{1,1}=[0.1 0.6 0.2]; +colors{2,1}=[0.4 0.1 0.4]; +colors{3,1}=[0.3 0.3 0.3]; +regcolors={}; +regcolors{1,1}=[0 0.1 0.4]; +regcolors{2,1}=[0 0.4 0.9]; +regcolors{3,1}=[0.5 0.5 1]; + +plotWindow=[-0.5 6]; +plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); +titles={'CS+ ','CS50 ','CS- '}; +sets={'set1','set2'}; + +%plot only in 2% +sel=sels{3}&~sels{5}; +cueResp=mean(cspActivity(sel,cueBins),2); +[~,sortOrder]=sort(cueResp); + +pn=0; +for cue=1:3 + for os=1:length(sets) + sn=sets{os}; + + activity = dffPSTH3{cue,os}(sel,:); + activity=activity(sortOrder,:); + pn=pn+1; + s=subplot(2,12,12+pn); + hold on; + + imagesc(binTimes(plotBins),[1 sum(sel)],activity(:,plotBins),[-2 2]);%[min(min(cspActivity)) max(max(cspActivity))]); + ylim([0.5 sum(sel)+0.5]) + plot([0 0],[0.5 sum(sel)+0.5],'color',colors{cue},'linewidth',0.75); + plot([2.5 2.5],[0.5 sum(sel)+0.5],':','color','k','linewidth',0.25); + set(gca,'ytick',[]); + + %colorbar; + title(titles{cue}); + if pn==1 + xlabel('seconds from odor onset'); + ylabel(sprintf('cue neurons in 2%% but not 5%% (%g)',sum(sel))); + end + colormap(s,map); + end +end +%plot in both 2% and 5% +sel=sels{3}&sels{5}; +cueResp=mean(cspActivity(sel,cueBins),2); +[~,sortOrder]=sort(cueResp); + +pn=0; +for cue=1:3 + for os=1:length(sets) + sn=sets{os}; + + activity = dffPSTH3{cue,os}(sel,:); + activity=activity(sortOrder,:); + pn=pn+1; + s=subplot(2,12,18+pn); + hold on; + + imagesc(binTimes(plotBins),[1 sum(sel)],activity(:,plotBins),[-2 2]);%[min(min(cspActivity)) max(max(cspActivity))]); + ylim([0.5 sum(sel)+0.5]) + plot([0 0],[0.5 sum(sel)+0.5],'color',colors{cue},'linewidth',0.75); + plot([2.5 2.5],[0.5 sum(sel)+0.5],':','color','k','linewidth',0.25); + set(gca,'ytick',[]); + + %colorbar; + title(titles{cue}); + if pn==1 + xlabel('seconds from odor onset'); + ylabel(sprintf('cue neurons in both (%g)',sum(sel))); + end + colormap(s,map); + end +end %% function A = makeA(discreteUnfixed,continuousPredictorValues) diff --git a/Code/ephysValue.m b/Code/ephysValue.m index 6d76acd..3b8284a 100644 --- a/Code/ephysValue.m +++ b/Code/ephysValue.m @@ -994,28 +994,10 @@ %% perform regression with reduced rank predictors preRun=true; -glmFile='valueGLM20220611'; +glmFile='valueGLM20230127'; tic - -%permute all possible combinations of odors for value shuffle -allodors=1:6; -perms6=NaN(90,6); -first2=nchoosek(allodors,2); -perm=0; -for f2=1:length(first2) - remaining=allodors(~ismember(allodors,first2(f2,:))); - next2=nchoosek(remaining,2); - for n2=1:length(next2) - last=remaining(~ismember(remaining,next2(n2,:))); - perm=perm+1; - perms6(perm,:)=[first2(f2,:) next2(n2,:) last]; - end -end -perms6all=perms6; - - if preRun load(glmFile); else @@ -1050,42 +1032,7 @@ As{sub} = At; rAs{sub} = At * bR(:,1:components); end - - %get original cue order for shuffle - firstBins={}; - cueOrder=[1 3 5 2 4 6]; - for cue=1:6 - firstBins{cue}=find(A(:,1+(cue-1)*binsPerCue)); - firstBins{cue}=[firstBins{cue} ones(length(firstBins{cue}),1)*cueOrder(cue)]; - end - alltrials=cat(1,firstBins{:}); - [sortedtrials,ordered]=sort(alltrials(:,1)); - originalOrder=alltrials(ordered,2); - - %make new shuffled predictor matrices - A3Sh={}; - for perm=1:length(perms6all)+1 - if perm>size(perms6all,1) %last model is with all cues weighted equally - newValue=ones(length(alltrials),1); - else - newOrder=NaN(length(alltrials),1); - for cue=1:6 - newOrder(originalOrder==cue)=perms6all(perm,cue); - end - - newValue=NaN(length(alltrials),1); - for cue=1:6 - originalValues=trialValues{session}(originalOrder==cue); - newValue(newOrder==cue)=mean(originalValues(randsample(sum(originalOrder==cue),sum(newOrder==cue),'true'))); - end - end - A3Sh{perm}=Xall3{session}; - for bin=1:binsPerCue - A3Sh{perm}(sortedtrials(:,1)-1+bin,bin)=newValue; - end - end - varExpLam=[]; for neuron=1:size(Yall{NS},2) NN=NN+1; @@ -1134,58 +1081,97 @@ varExp(NN,1+sub) = 1- var(y-pred)/var(y); end - - %value instead of cue identitiy - A3 = Xall3{NS}; - rA3 = A3 * bR3(:,1:components); - pred=NaN(size(y,1),1); - for fold=1:folds - train=trains{fold}; - test=train==0; - fitK=lassoglm(rA3(train,:),y(train),'normal','alpha',0.5,'lambda',thisLambda); - kernel=bR3(:,1:components)*fitK; - pred(test)=A3(test,:)*kernel; - end - varExp(NN,3+sub) = 1- var(y-pred)/var(y); - - for perm=1:size(perms6all,1)+1 - - rA3 = A3Sh{perm} * bR3(:,1:components); - trains = getfolds(rA3,folds,binsPerTrial); - - % if ismember(perm,topModels) - % display(perms6all(perm,:));display([combinedValues(1:10) newValue(1:10)]); - % end - - - pred=NaN(size(y,1),1); - for fold=1:folds - train=trains{fold}; - test=train==0; - fitK=lassoglm(rA3(train,:),y(train),'normal','alpha',0.5,'lambda',thisLambda); - kernel=bR3(:,1:components)*fitK; - pred(test)=A3Sh{perm}(test,:)*kernel; - end - varExpValueSh(NN,3+perm) = 1- var(y-pred)/var(y); - end - + end fprintf('Session #%d \n',NS); end -varExpValueSh(:,1)=varExp(:,1); -varExpValueSh(:,2)=varExp(:,2); -varExpValueSh(:,3)=varExp(:,5); +save('newGLMFile.mat','varExp','kernels','lambdaVal','predF','-v7.3'); + +end +toc + + +%% perform value analysis with a ton of shuffles +runglm=false; +%generate all value shuffles +valueShuffles = []; +valueVersion = []; +valueLevels = [1 1 0.7 0.7 0.1 0.1;1 1 0.5 0.5 0 0;1 1 1 1 0 0;1 1 1 0 0 0;1 1 0 0 0 0;1 0 0 0 0 0;1 1 1 1 1 0;1 1 1 1 1 1]; +for vl=1:size(valueLevels,1) +allCombos = perms(valueLevels(vl,:)); +allCombos = unique(allCombos,'rows'); +valueShuffles = cat(1,valueShuffles,allCombos); +valueVersion = cat(1,valueVersion,vl*ones(size(allCombos,1),1)); +end +valueShuffles=fliplr(valueShuffles); + +folds=4; +NNstart=0; +%varExpValueNew=NaN(totalNeurons,length(valueShuffles)); +%kernels=NaN(totalNeurons,31); +binsPerCue=diff(windows{1})/binSize; + +opts.alpha=0.5; +glmopts=glmnetSet(opts); -save('newGLMFile.mat','varExp','varExpValueSh','kernels','lambdaVal','predF','-v7.3'); +if runglm +sessInd=find(includedSessions); +tic +for NS=1:sum(includedSessions) + session=sessInd(incses); + sessionFolder=fullfile(direc,sessionSubject{session},sessionDate{session}); + cueTimes = readNPY(fullfile(sessionFolder,'cues.times.npy')); + rewprob = readNPY(fullfile(sessionFolder,'cues.rewardProbability.npy')); + odorset = readNPY(fullfile(sessionFolder,'cues.odorSet.npy')); + cueOrder=zeros(length(cueTimes),1); + cueOrder(rewprob==1&odorset==1)=1; + cueOrder(rewprob==1&odorset==2)=2; + cueOrder(rewprob==0.5&odorset==1)=3; + cueOrder(rewprob==0.5&odorset==2)=4; + cueOrder(rewprob==0&odorset==1)=5; + cueOrder(rewprob==0&odorset==2)=6; + + A = Xall{NS}(:,[1:binsPerCue size(Xall{NS},2)-5:size(Xall{NS},2)]); + trains = getfolds(A,folds,binsPerTrial); + for s=length(valueShuffles)-6:length(valueShuffles) %1:length(valueShuffles) + for c=1:length(cueTimes) + cueValue=valueShuffles(s,cueOrder(c)); + dm=cueValue*eye(binsPerCue); + A(c*binsPerTrial-binsPerCue+1:c*binsPerTrial,1:binsPerCue)=dm; + end + NN=NNstart; + for neuron=1:size(Yall{NS},2) + NN=NN+1; + y=Yall{NS}(:,neuron); + thisLambda=lambdaVal(NN,1); + %cross-validated variance explained + pred=NaN(size(y,1),1); + for fold=1:folds + train=trains{fold}; + test=train==0; + %kernel=lassoglm(A(train,:),y(train),'normal','alpha',0.5,'lambda',thisLambda); + fit = glmnet(A(train,:),y(train),[],glmopts); + prediction = glmnetPredict(fit,A(test,:)); + pred(test)=prediction(:,end); + end + varExpValueNew(NN,s) = 1- var(y-pred)/var(y); + kernel = fit.beta(:,end); + if s==1 kernels(NN,:)=kernel; end + end + end + NNstart=NN; + fprintf('Session #%d \n',NS); end toc +save('valueGLMFile.mat','varExp','varExpValueNew','kernels','lambdaVal','predF','-v7.3'); +end -%% value coding analysis 90 with untuned +%% value coding analysis improvCutoff=0.02; overallCutoff=0.02; improvement=varExp(:,1)-varExp; @@ -1196,114 +1182,150 @@ cueK = improvement(:,2)>improvCutoff; %if best model is at least cutoff better than no cue kernels totalCueVariance=varExp(:,1)-varExp(:,2); %best model compared to no cue model -valueImp=varExp(:,1)-varExpValueSh; %how much better the best model is than value, shuffles, and one kernel +valueImp=varExp(:,1)-varExpValueNew; %how much better the best model is than value, shuffles, and one kernel %correct it so you can't be worse than no cues at all (should I do this?) for n=1:length(valueImp) valueImp(n,valueImp(n,:)>totalCueVariance(n))=totalCueVariance(n); end -%value specific variance -%whatever is accounted for by value kernel -cueVariance=totalCueVariance - valueImp(:,[4:end]); -[~,bestModel]=max(cueVariance,[],2); -bestModel(isnan(cueVariance(:,1)))=NaN; -%region -neuronRegionOlf=neuronRegionAdj; -neuronRegionOlf(ismember(neuronRegionAdj,{'EPd','PIR'}))={'OLF'}; -%neuronRegionOlf(ismember(neuronRegionAdj,{'EPd'}))={'OLF'}; - -division=neuronRegionOlf; -%division=neuronSubregion; - -regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; -%regions={'MOs','ACAd','FRP','PL','ILA','ORBl','ORBvl','ORBm','DP','TTd','AON'}; -vcolors={[0 0.7 1],[0.6 0.6 0.6]}; +includedModels = [91:length(valueShuffles)]; +%includedModels = [1 91 2:90 92:length(valueShuffles)]; +%includedModels = [1:90 181:length(valueShuffles)]; +%includedModels = [1:length(valueShuffles)]; +shuffleLabels = cell(length(includedModels),1); +for sl=1:length(includedModels) + m=includedModels(sl); + shuffleLabels{sl} = sprintf('%g,%g,%g,%g,%g,%g',valueShuffles(m,1),valueShuffles(m,2),valueShuffles(m,3),valueShuffles(m,4),valueShuffles(m,5),valueShuffles(m,6)); +end +[~,bestModel]=max(varExpValueNew(:,includedModels),[],2); +% [~,bestModel]=max(varExpValueNew,[],2); +% bestModel(bestModel>90)=bestModel(bestModel>90)-90; +totalModels = length(includedModels); figure; -cueNames={'+','+','50','50','-','-'}; -odorValues=[0.5 0.5 0.37 0.37 0.05 0.05]; -vnames={}; -shuffVals=NaN(90,6); -for n=1:length(perms6all) - vnames{n,1}=append(cueNames{perms6all(n,:)==1},'/',cueNames{perms6all(n,:)==2},', ',cueNames{perms6all(n,:)==3},'/',cueNames{perms6all(n,:)==4},', ',cueNames{perms6all(n,:)==5},'/',cueNames{perms6all(n,:)==6}); - for odor=1:6 - shuffVals(n,odor)=odorValues(perms6all(n,:)==odor); - end -end +corrWVal=corrcoef(valueShuffles(includedModels,:)'); +corrWVal(isnan(corrWVal))=0; +[valCorr,valRank] = sort(corrWVal(:,1),'descend'); +[~,mdlPosition] = sort(valRank); -subplot(3,1,1); +subplot(5,1,1); hold on -frequency=histcounts(bestModel(cueK&~behavior&~isnan(bestModel)),0.5:1:91.5,'normalization','probability'); -outliers=frequency>0.017; -topModels=find(outliers); -b=bar([1:max(bestModel)-1 max(bestModel)+1],frequency,'facecolor',[0 0.2 0.4],'linewidth',0.01); +frequency=histcounts(bestModel(cueK&~behavior&~isnan(bestModel)),0.5:1:totalModels+0.5,'normalization','probability'); +%b=bar(1:max(bestModel),frequency,'facecolor',[0 0.2 0.4],'linewidth',0.01); +b=bar(mdlPosition,frequency,'facecolor',[0 0.2 0.4],'linewidth',0.01); + b.FaceColor='flat'; -for e=1:length(topModels) -b.CData(topModels(e),:)=[0.7 0 1]; %trial type +for e=2:17 +%b.CData(valRank(e),:)=[0.7 0 1]; %trial type +b.CData(e,:)=[0.7 0 1]; %trial type end b.CData(1,:)=[0 0.7 1]; %value -b.CData(end,:)=[0.6 0.6 0.6]; %non-specific - +b.CData(mdlPosition(end),:)=[0.6 0.6 0.6]; %non-specific +%b.CData(end,:)=[0.6 0.6 0.6]; %non-specific +ylabel('fraction of cue cells'); +chance=1/totalModels; +plot([0 totalModels+1],[chance chance],':','color','k','linewidth',0.5); +xlim([0 totalModels+1]); +ylim([0 0.2]); +%xticks(mdlPosition(topModels(1:end-1))); +xticks([]); +yticks(0:0.1:0.3); +xlabel('cue coding models, sorted by similarity to value model'); +subplot(4,5,6); +hold on +scatter(valCorr(2:17,1),frequency(valRank(2:17)),24,[0.7 0 1],'filled'); +scatter(valCorr(1,1),frequency(valRank(1)),24,[0 0.7 1],'filled'); +scatter(valCorr(18:end,1),frequency(valRank(18:end)),24,[0 0.2 0.4],'filled'); +scatter(valCorr(mdlPosition(end),1),frequency(end),24,[0.6 0.6 0.6],'filled'); +plot([-1 1],[1/length(includedModels) 1/length(includedModels)],':','color','k','linewidth',0.5); +ylim([0 0.2]); ylabel('fraction of cue cells'); +xlabel('correlation with value model'); +yticks(0:0.1:0.3); -chance=1/90 * sum(bestModel(cueK&~behavior&~isnan(bestModel))<91)/sum(cueK&~behavior&~isnan(bestModel)); -plot([0 91],[chance chance],':','color','k','linewidth',0.5); -xlim([0 93]); -ylim([0 0.35]); -xticks(topModels(1:end-1)); +%trial type +tt=valueShuffles(includedModels,1)==valueShuffles(includedModels,2) &... + valueShuffles(includedModels,3)==valueShuffles(includedModels,4) &... + valueShuffles(includedModels,5)==valueShuffles(includedModels,6); +subplot(4,5,7); +hold on +scatter(corrWVal(tt==0,1),frequency(tt==0),24,[0 0 0],'filled'); +scatter(corrWVal(tt,1),frequency(tt),24,[1 0 0.3],'filled'); +plot([-1 1],[1/length(includedModels) 1/length(includedModels)],':','color','k','linewidth',0.5); +ylim([0 0.2]); +ylabel('fraction of cue cells'); +xlabel('correlation with value model'); yticks(0:0.1:0.3); -%xtickangle(45); -%xticklabels(vnames(outliers)); + category=9*ones(length(improvement),1); -category(cueK&~behavior)=8; %cue neurons -%categorize neurons by their best model if it's in the top models -for bm=1:length(topModels) -category(cueK&~behavior&bestModel==topModels(bm))=bm; %cue, odor -end +category(cueK&~behavior)=4; %cue neurons, other odor coding schemes +category(cueK&~behavior&bestModel==1)=1; %value +category(cueK&~behavior&ismember(bestModel,valRank(2:17)))=2; %value like +category(cueK&~behavior&bestModel==max(bestModel))=3; %untuned + + + + +subplot(4,5,8); +hold on +valcand=ismember(bestModel,[1 2]); +scatter(varExpValueNew(valcand,1),varExpValueNew(valcand,91)); +plot([0 1],[0 1],'color','k'); +xlabel('variance explained with 1,0.7,0.1'); +ylabel('variance explained with 1,0.5,0'); +title('value cells (selected with either)'); % visual depiction of shuffles colormap('summer'); -subplot(5,1,3); -imagesc(1:45,1:6,shuffVals(1:45,:)',[0 0.5]); -xticks(topModels); +subplot(5,1,4); +imagesc(1:length(includedModels),1:6,valueShuffles(91:end,:)',[0 1]); +xticks([]); yticks(1:6); yticklabels({'CS+','CS+','CS50','CS50','CS-','CS-'}); title('Value assigned to each odor in each shuffle'); +colorbar; -subplot(5,1,4); -imagesc(46:90,1:6,shuffVals(46:end,:)',[0 0.5]); -xticks(topModels); +% visual depiction of shuffles, reordered by correlation +subplot(5,1,5); +includedShuffles = valueShuffles(91:end,:)'; + +imagesc(1:length(includedModels),1:6,includedShuffles(:,valRank),[0 1]); +xticks([]); yticks(1:6); yticklabels({'CS+','CS+','CS50','CS50','CS-','CS-'}); - - -subplot(5,2,9); -imagesc(1:45,1:6,shuffVals(1:45,:)',[0 0.5]); -xticks(topModels); +title('Value assigned to each odor in each shuffle, ordered by similarity to value'); colorbar; -%% get proportions of value and trial type and compare between regions +%% get proportions of value and meaning and compare between regions regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; -reggrps=[1:3]; +reggrps=[2:4]; +neuronGroupAdj = neuronGroup + 1; +neuronGroupAdj(neuronGroupAdj==5)=1; +pcutoff=0.05/36; % regions={'ACA','FRP','PL','ILA','ORB','CP','ACB'}; % reggrps=[2 4]; - +neuronRegionOlf=neuronRegionAdj; +neuronRegionOlf(ismember(neuronRegionAdj,{'EPd','PIR'}))={'OLF'}; allReg=unique(neuronRegionOlf); +order=[9 10 5 2 3 1 7 12 8 11 6 13 4]; +allReg=allReg(order); + regInd=NaN(length(regions),1); for reg=1:length(regions) regInd(reg,1)=find(ismember(allReg,regions(reg))); end + %fit all together, allowing random effect of session to have more power catcolors={[0 0.7 1],[0.7 0 1],[0.6 0.6 0.6],[0.05 0.25 0.45]}; -titles={'value','trial type','non-specific','other selective'}; -cats={1,2:6,7,8}; -for ct=1:4 +titles={'value','value-like','untuned','other'}; +cats={1,2,3,4}; +for ct=1:length(cats) %region barData=zeros(2,5); subplot(4,4,ct); @@ -1318,13 +1340,43 @@ barData(2,reg)=sum(othersel®sel)/sum(regsel); end tbl = table(categorical(neuronSession(regionnumber>0)), categorical(regionnumber(regionnumber>0)), sel(regionnumber>0), 'VariableNames', {'session', 'region','val'}); - lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial'); - [ypred,ypredCI]=predict(lm,'conditional',false); + lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial','dummyvarcoding','reference'); + + [ypred,ypredCI,df]=predict(lm,'conditional',false); + + Hall=[]; + Htemp = [0 0 0 0 0 0 0 0 0 0 0 0 0]; + for reg1 = 1:length(regions)-1 + for reg2 = reg1+1:length(regions) + H=Htemp; + H(reg1+4)=1; + H(reg2+4)=-1; + Hall=cat(1,Hall,H); + end + end + greater=zeros(length(regions)); + for reg1 = 1:length(regions) + for reg2 = 1:length(regions) + H=Htemp; + H(reg1+4)=1; + H(reg2+4)=-1; + greater(reg1,reg2)=coefTest(lm,H); + end + end + %Hall=unique(Hall,'rows'); + pVals=NaN(length(Hall),1); + for p=1:length(pVals) + pVals(p,1) = coefTest(lm,Hall(p,:)); + end barData=barData(:,regInd); b=bar(barData','stacked'); - b(2).FaceColor=[0.85 0.85 0.85]; - b(1).FaceColor=catcolors{ct}; + %b(3).FaceColor=[0.4 0 0.5]; %odor + b(2).FaceColor=[0.85 0.85 0.85]; %meaning + b(1).FaceColor=catcolors{ct}; %value + %b(4).FaceColor=[0 0 0.2]; %odor + % b(5).FaceColor=[0.9 0.2 0]; + %b(3).FaceColor=[1 1 1]; upperCI=[]; lowerCI=[]; estmean=[]; @@ -1354,16 +1406,16 @@ customMap=cat(1,[1 1 1],customMap); %by region - greater=zeros(length(regions)); + %greater=zeros(length(regions)); frac=zeros(length(regions)); for r1=1:length(regions) for r2=1:length(regions) frac(r1,r2)=barData(1,r1)-barData(1,r2); - greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); + %greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end s=subplot(4,4,4+ct); - imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); + imagesc(frac .* (greater0), [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); cb=colorbar; @@ -1387,24 +1439,39 @@ othersel=ismember(category,1:8)&~sel; regionnumber=zeros(totalNeurons,1); for reg=1:4 - regsel=ismember(neuronGroup,reg); + regsel=ismember(neuronGroupAdj,reg); regionnumber(regsel)=reg; barData(1,reg)=sum(sel®sel)/sum(regsel); barData(2,reg)=sum(othersel®sel)/sum(regsel); end tbl = table(categorical(neuronSession(regionnumber>0)), categorical(regionnumber(regionnumber>0)), sel(regionnumber>0), 'VariableNames', {'session', 'region','val'}); - lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial'); + lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial','dummyvarcoding','reference'); [ypred,ypredCI]=predict(lm,'conditional',false); - + + Htemp = [0 0 0 0]; + greater=zeros(length(reggrps)); + for reg1 = 1:length(reggrps) + for reg2 = 1:length(reggrps) + H=Htemp; + H(reg1+1)=1; + H(reg2+1)=-1; + greater(reg1,reg2)=coefTest(lm,H); + end + end + barData=barData(:,reggrps); b=bar(barData','stacked'); + %b(3).FaceColor=[0.4 0 0.5]; %odor b(2).FaceColor=[0.85 0.85 0.85]; - b(1).FaceColor=catcolors{ct}; + b(1).FaceColor=catcolors{ct}; %value + %b(4).FaceColor=[0 0 0.2]; %odor + % b(5).FaceColor=[0.9 0.2 0]; + %b(3).FaceColor=[1 1 1]; upperCI=[]; lowerCI=[]; estmean=[]; for reg=1:length(reggrps) - regsel=find(ismember(neuronGroup,reggrps(reg)),1); + regsel=find(ismember(neuronGroupAdj,reggrps(reg)),1); estmean(1,reg)=ypred(regsel); upperCI(1,reg)=ypredCI(regsel,2); lowerCI(1,reg)=ypredCI(regsel,1); @@ -1412,6 +1479,19 @@ errorbar(1:length(reggrps),estmean(1,:),estmean(1,:)-lowerCI(1,:),upperCI(1,:)-estmean(1,:),'o','linewidth',0.75,'color','k'); + + Htemp = [0 0 0 0]; + greater=zeros(length(reggrps)); + for reg1 = 1:length(reggrps) + for reg2 = 1:length(reggrps) + H=Htemp; + H(reg1+1)=1; + H(reg2+1)=-1; + greater(reg1,reg2)=coefTest(lm,H); + end + end + + xlim([0.5 7+0.5]); %legend('Cues','Cues,Licks','Licks','Licks,Reward','Reward','Cues,Reward','Cues,Licks,Reward'); xticks(1:length(reggrps)); @@ -1429,25 +1509,25 @@ customMap=cat(1,[1 1 1],customMap); %by region - greater=zeros(length(reggrps)); + %greater=zeros(length(reggrps)); frac=zeros(length(reggrps)); for r1=1:length(reggrps) for r2=1:length(reggrps) frac(r1,r2)=barData(1,r1)-barData(1,r2); - greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); + %greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end s=subplot(4,4,12+ct); - imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); + imagesc(frac .* (greater<0.05/3 & frac>0), [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); cb=colorbar; set(cb,'ticks',[0:0.1:0.3]); xticks(1:length(reggrps)); - xticklabels(regionGroupNames(reggrps)); + xticklabels(regionGroupNames(reggrps-1)); xtickangle(45); yticks(1:length(reggrps)); - yticklabels(regionGroupNames(reggrps)); + yticklabels(regionGroupNames(reggrps-1)); title(titles{ct}); if ct==1 ylabel('increase in this region...'); @@ -1457,14 +1537,19 @@ end -%% get proportions of value and trial type out of cue cells +%% get proportions of value and meaning out of cue cells regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; -reggrps=[1:3]; +reggrps=[2:4]; +neuronGroupAdj = neuronGroup + 1; +neuronGroupAdj(neuronGroupAdj==5)=1; +pcutoff=0.05/36; %regions={'ACA','FRP','PL','ILA','ORB','CP','ACB'}; %reggrps=[2 4]; allReg=unique(neuronRegionOlf); +order=[9 10 5 2 3 1 7 12 8 11 6 13 4]; +allReg=allReg(order); regInd=NaN(length(regions),1); for reg=1:length(regions) regInd(reg,1)=find(ismember(allReg,regions(reg))); @@ -1472,9 +1557,9 @@ %fit all together, allowing random effect of session to have more power catcolors={[0 0.7 1],[0.7 0 1],[0.6 0.6 0.6],[0.05 0.25 0.45]}; -titles={'value','trial type','untuned','other selective'}; -cats={1,2:6,7,8}; -for ct=1:4 +titles={'value','value-like','untuned','other'}; +cats={1,2,3,4}; +for ct=1:length(cats) %region barData=zeros(2,5); subplot(4,4,ct); @@ -1492,10 +1577,42 @@ lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial'); [ypred,ypredCI]=predict(lm,'conditional',false); + + Hall=[]; + Htemp = [0 0 0 0 0 0 0 0 0 0 0 0 0]; + for reg1 = 1:length(regions)-1 + for reg2 = reg1+1:length(regions) + H=Htemp; + H(reg1+4)=1; + H(reg2+4)=-1; + Hall=cat(1,Hall,H); + end + end + greater=zeros(length(regions)); + for reg1 = 1:length(regions) + for reg2 = 1:length(regions) + H=Htemp; + H(reg1+4)=1; + H(reg2+4)=-1; + greater(reg1,reg2)=coefTest(lm,H); + end + end + + %Hall=unique(Hall,'rows'); + pVals=NaN(length(Hall),1); + for p=1:length(pVals) + pVals(p,1) = coefTest(lm,Hall(p,:)); + end + + barData=barData(:,regInd); b=bar(barData','stacked'); - b(2).FaceColor=[0.85 0.85 0.85]; - b(1).FaceColor=catcolors{ct}; + %b(3).FaceColor=[0.4 0 0.5]; %odor + b(2).FaceColor=[0.85 0.85 0.85]; %meaning + b(1).FaceColor=catcolors{ct}; %value + %b(4).FaceColor=[0 0 0.2]; %odor + % b(5).FaceColor=[0.9 0.2 0]; + %b(3).FaceColor=[1 1 1]; upperCI=[]; lowerCI=[]; estmean=[]; @@ -1505,7 +1622,7 @@ upperCI(1,reg)=ypredCI(regsel,2); lowerCI(1,reg)=ypredCI(regsel,1); end - + errorbar(1:length(regions),estmean(1,:),estmean(1,:)-lowerCI(1,:),upperCI(1,:)-estmean(1,:),'o','linewidth',0.75,'color','k'); xlim([0.5 length(regions)+0.5]); @@ -1525,16 +1642,16 @@ customMap=cat(1,[1 1 1],customMap); %by region - greater=zeros(length(regions)); + %greater=zeros(length(regions)); frac=zeros(length(regions)); for r1=1:length(regions) for r2=1:length(regions) frac(r1,r2)=barData(1,r1)-barData(1,r2); - greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); + %greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end s=subplot(4,4,4+ct); - imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); + imagesc(frac .* (greater0), [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); cb=colorbar; @@ -1549,6 +1666,7 @@ ylabel('increase in this region...'); xlabel('over this region'); end + %region group barData=zeros(2,3); @@ -1558,7 +1676,7 @@ othersel=ismember(category,1:8)&~sel; regionnumber=zeros(totalNeurons,1); for reg=1:4 - regsel=ismember(neuronGroup,reg)&category<9; + regsel=ismember(neuronGroupAdj,reg)&category<9; regionnumber(regsel)=reg; barData(1,reg)=sum(sel®sel)/sum(regsel); barData(2,reg)=sum(othersel®sel)/sum(regsel); @@ -1568,15 +1686,32 @@ lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial'); [ypred,ypredCI]=predict(lm,'conditional',false); + + Htemp = [0 0 0 0]; + greater=zeros(length(reggrps)); + for reg1 = 1:length(reggrps) + for reg2 = 1:length(reggrps) + H=Htemp; + H(reg1+1)=1; + H(reg2+1)=-1; + greater(reg1,reg2)=coefTest(lm,H); + end + end + + barData=barData(:,reggrps); b=bar(barData','stacked'); + %b(3).FaceColor=[0.4 0 0.5]; %odor b(2).FaceColor=[0.85 0.85 0.85]; - b(1).FaceColor=catcolors{ct}; + b(1).FaceColor=catcolors{ct}; %value + %b(4).FaceColor=[0 0 0.2]; %odor + % b(5).FaceColor=[0.9 0.2 0]; + %b(3).FaceColor=[1 1 1]; upperCI=[]; lowerCI=[]; estmean=[]; for reg=1:length(reggrps) - regsel=find(ismember(neuronGroup(category<9),reggrps(reg)),1); + regsel=find(ismember(neuronGroupAdj(category<9),reggrps(reg)),1); estmean(1,reg)=ypred(regsel); upperCI(1,reg)=ypredCI(regsel,2); lowerCI(1,reg)=ypredCI(regsel,1); @@ -1587,7 +1722,7 @@ xlim([0.5 7+0.5]); %legend('Cues','Cues,Licks','Licks','Licks,Reward','Reward','Cues,Reward','Cues,Licks,Reward'); xticks(1:length(reggrps)); - xticklabels(regionGroupNames(reggrps)); + xticklabels(regionGroupNames(reggrps-1)); xtickangle(45); ylim([0 1]); @@ -1601,25 +1736,25 @@ customMap=cat(1,[1 1 1],customMap); %by region - greater=zeros(length(reggrps)); + %greater=zeros(length(reggrps)); frac=zeros(length(reggrps)); for r1=1:length(reggrps) for r2=1:length(reggrps) frac(r1,r2)=barData(1,r1)-barData(1,r2); - greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); + %greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end s=subplot(4,4,12+ct); - imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); + imagesc(frac .* (greater<0.05/3 & frac>0), [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); cb=colorbar; set(cb,'ticks',[0:0.1:0.3]); xticks(1:length(reggrps)); - xticklabels(regionGroupNames(reggrps)); + xticklabels(regionGroupNames(reggrps-1)); xtickangle(45); yticks(1:length(reggrps)); - yticklabels(regionGroupNames(reggrps)); + yticklabels(regionGroupNames(reggrps-1)); title(titles{ct}); if ct==1 ylabel('increase in this region...'); @@ -1630,8 +1765,8 @@ %% value models schematic -examplePerms=[1 33]; -exampleNeurons=[9]; +examplePerms=[91 123]; +exampleNeurons=[3 10 20]; plotWindow=[0 2.5]; plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); @@ -1645,13 +1780,6 @@ colors{3,1}=[0.3 0.3 0.3]; specs={'-','--'}; - -%smoothing filter for licks -smoothbinsLick=25; %number of previous bins used to smooth -halfnormal=makedist('HalfNormal','mu',0,'sigma',20); -filterweightsLick=pdf(halfnormal,0:smoothbinsLick); - -components=10; folds=4; neuronNumber=[1:length(neuronSession)]'; @@ -1686,268 +1814,106 @@ ylabel('z-score'); end title(num2str(exampleNeurons(n))); - + %get model fits session=neuronSession(nn); sessionNeuron=find(neuronNumber(neuronSession==session)==nn); sessionFolder=fullfile(direc,sessionSubject{session},sessionDate{session}); - lickTimes = readNPY(fullfile(sessionFolder,'licks.times.npy')); - reward = readNPY(fullfile(sessionFolder,'rewards.times.npy')); cueTimes = readNPY(fullfile(sessionFolder,'cues.times.npy')); rewprob = readNPY(fullfile(sessionFolder,'cues.rewardProbability.npy')); - odorset = readNPY(fullfile(sessionFolder,'cues.odorSet.npy')); - csp1 = cueTimes(rewprob==1&odorset==1); - csp2 = cueTimes(rewprob==1&odorset==2); - csf1 = cueTimes(rewprob==0.5&odorset==1); - csf2 = cueTimes(rewprob==0.5&odorset==2); - csm1 = cueTimes(rewprob==0&odorset==1); - csm2 = cueTimes(rewprob==0&odorset==2); - - %first lick in bout - ibi=0.5; - boutStart = lickTimes(cat(1,1,1+find(diff(lickTimes)>ibi))); - - binEdgesSession=cueTimes(1)-5:binSizeSession:cueTimes(end)+10; - binTimesSession=cueTimes(1)-5+binSizeSession/2:binSizeSession:cueTimes(end)+10-binSizeSession/2; - binnedLicksGLM=histcounts(lickTimes,binEdgesSession) / binSizeSession; - lickRateGLM=NaN(length(binnedLicksGLM),1); - for l=1:length(binnedLicksGLM) - lickRateGLM(l,1)=sum(binnedLicksGLM(l-min([l-1 smoothbinsLick]):l).*fliplr(filterweightsLick(1:min([l smoothbinsLick+1]))))/sum(filterweightsLick(1:min([l smoothbinsLick+1]))); - end - - binnedLicksGLM=NaN(length(cueTimes),length(binEdges)-1); - binTimesGLM=NaN(length(cueTimes),length(binEdges)-1); - for trial = 1:length(cueTimes) - binTimesRel = (binTimesSession - cueTimes(trial)); - for bin=1:length(binTimes) - binnedLicksGLM(trial,bin)=mean(lickRateGLM(binTimesRel>=binEdges(bin)&binTimesRel1 + xticks([]); + end + if n==1 + title(sprintf('shuffle %d',examplePerms(s))); end + end - binnedLicksGLM(isnan(binnedLicksGLM))=0; %this only happens if bins go beyond end of session, so extremely rare + +end - vector1={}; - includedBins=find(binTimes>analysisWindow(1) & binTimes1 - xticks([]); - end - if n==1 - title(sprintf('shuffle %d',examplePerms(permn))); - end - - end - -end - - -%% activity of cue cell types and correlations -figure; -map=redblue(256); -CL=[0 1]; cueWindow=[0 2.5]; cueBins=binTimes>=cueWindow(1) & binTimes<=cueWindow(2); -cueCorrMatrix=NaN(12,12,totalNeurons); -cueCorrSets=NaN(6,6,totalNeurons); -setCorr=NaN(totalNeurons,3); -orders=[1 2 3;2 3 1;3 1 2]; - -for NN=1:length(neuronRegionOlf) - cueVector=NaN(sum(cueBins),12); - cue3Vector=NaN(sum(cueBins)*3,6); - v=0; - u=0; - for os=1:2 - for condition=1:6 - v=v+1; - cueVector(:,v)=PSTH6{condition,os}(NN,cueBins); - end - - end - %cueVector=cueVector-mean(cueVector,2); - corrmat=corrcoef(cueVector); - - cueCorrMatrix(:,:,NN)=corrmat; -end - - %heatmaps plotWindow=[-0.5 2.5]; plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); -sortWindow=[0 1.5]; +sortWindow=[0 2.5]; sortBins=binTimes>=sortWindow(1) & binTimes<=sortWindow(2); - - colors{1,1}=[0.1 0.6 0.2]; colors{2,1}=[0.4 0.1 0.4]; colors{3,1}=[0.3 0.3 0.3]; - -for ct=1 - sel=category==ct; -thisActivity=PSTH3{1,1}(sel,:); -cueResp=mean(thisActivity(:,sortBins),2); -[~,sortOrder]=sort(cueResp); -pn=0; -activity=[]; -for cue=1:3 - for os=1:2 - - pn=pn+1; - - ax(1)=subplot(2,20,20+(ct-1)*6+pn); - colormap(ax(1),map); - hold on; - - - activity = PSTH3{cue,os}(sel,:); - activity = activity(sortOrder,:); - - - - imagesc(binTimes(plotBins),[1 length(activity)],activity(:,plotBins),[-3 3]);%[min(min(cspActivity)) max(max(cspActivity))]); - ylim([0.5 size(activity,1)+0.5]) - plot([0 0],[0.5 sum(sel)+0.5],'color',colors{cue},'linewidth',0.75); - set(gca,'ytick',[]); - if pn==1 xlabel('seconds from cue'); end - if pn>1 xticks([]); end - - - - end - - - -end -end - - - pn=0; activity=[]; -sel=ismember(category,2:6); -thisActivity=PSTH3{1,1}(sel,:); +thisActivity=PSTH3{1,1}; cueResp=mean(thisActivity(:,sortBins),2); - -regionCodeOpp=8-category; -sortcrit=[regionCodeOpp(sel) cueResp]; +cats={1,2,3}; + +% regions={'ALM','MOs','ACA','FRP','PL','ILA','ORB','CP','ACB','DP','TTd','AON','OLF'}; +% +% neuronRegionOlf=neuronRegionAdj; +% neuronRegionOlf(ismember(neuronRegionAdj,{'EPd','PIR'}))={'OLF'}; +% regionCode=NaN(length(neuronRegionOlf),1); +% for NN=1:length(neuronRegionOlf) +% regionCode(NN,1)=find(ismember(regions,neuronRegionOlf(NN))); +% end +% regionCodeOpp=12-regionCode; +regionCodeOpp=ones(totalNeurons,1); + +for ct=1:length(cats) +sel=ismember(category,cats{ct}); +sortcrit=[regionCodeOpp(sel) cueResp(sel)]; [~,sortOrder]=sortrows(sortcrit,[1 2]); +pn=0; for cue=1:3 for os=1:2 pn=pn+1; - ax(1)=subplot(2,20,27+(ct-1)*6+pn); + ax(1)=subplot(3,15,(ct-1)*15+pn); colormap(ax(1),map); hold on; @@ -1961,6 +1927,7 @@ ylim([0.5 size(activity,1)+0.5]) plot([0 0],[0.5 sum(sel)+0.5],'color',colors{cue},'linewidth',0.75); set(gca,'ytick',[]); + if pn==1 ylabel(sprintf('n=%d',sum(sel))); end if pn==1 xlabel('seconds from cue'); end if pn>1 xticks([]); end @@ -1972,96 +1939,59 @@ end -pn=0; -activity=[]; -sel=ismember(category,7:8); -thisActivity=PSTH3{1,1}(sel,:); -cueResp=mean(thisActivity(:,sortBins),2); +colors{1,1}=[0.1 0.6 0.2]; +colors{2,1}=[0.4 0.1 0.4]; +colors{3,1}=[0.3 0.3 0.3]; +directions=[1 -1]; +specs={'-','--'}; +diractivity = mean(cat(3,PSTH6{(1-1)*2+2,1},PSTH6{(1-1)*2+1,2}),3); +direction=sign(max(diractivity,[],2)-abs(min(diractivity,[],2))); -regionCodeOpp=8-category; -sortcrit=[regionCodeOpp(sel) cueResp]; -[~,sortOrder]=sortrows(sortcrit,[1 2]); -for cue=1:3 - for os=1:2 - - pn=pn+1; + for d=1:2 + subplot(6,3,3+(ct-1)*6+(d-1)*3); - ax(1)=subplot(2,20,34+(ct-1)*6+pn); - colormap(ax(1),map); hold on; + for cue=1:3 + for os=1:2 + + activity = PSTH6{(cue-1)*2+os,os}; + + %get values + psth=nanmean(activity(sel&direction==directions(d),plotBins)); + sem=nanste(activity(sel&direction==directions(d),plotBins),1); %calculate standard error of the mean + up=psth+sem; + down=psth-sem; + + %plotting + plot(binTimes(plotBins),psth,specs{os},'Color',colors{cue,1},'linewidth',0.75); + %patch([xvals,xvals(end:-1:1)],[up,down(end:-1:1)],regcolors{reg,1},'EdgeColor','none');alpha(0.2); + + + end + + end + if d==1 text(0.1,1.5,num2str(sum(sel&direction==directions(d)))); end + if d==2 text(0.1,-0.45,num2str(sum(sel&direction==directions(d)))); end + if reg==1 & d==2 + xlabel('seconds from odor onset'); + ylabel('z-score'); + else + xticks([]); + end - activity = PSTH3{cue,os}(sel,:); - activity = activity(sortOrder,:); - - - - imagesc(binTimes(plotBins),[1 length(activity)],activity(:,plotBins),[-3 3]);%[min(min(cspActivity)) max(max(cspActivity))]); - ylim([0.5 size(activity,1)+0.5]) - plot([0 0],[0.5 sum(sel)+0.5],'color',colors{cue},'linewidth',0.75); - set(gca,'ytick',[]); - if pn==1 xlabel('seconds from cue'); end - if pn>1 xticks([]); end - - - + if reg>1 yticks([]); end + plot([0 0],[-1 2.5],'color','k','linewidth',0.75); + plot(plotWindow,[0 0],':','color','k','linewidth',0.5); + if d==1 axis([plotWindow -0.25 2.5]); end + if d==2 axis([plotWindow -1 0.2]); end end - - - -end -titles={vnames{outliers(1:end-1)},'non-specific','odor'}; -for ct=1:8 -ax(2)=subplot(4,8,8+ct); -sel=category==ct; -imagesc(nanmean(cueCorrMatrix(:,:,sel),3),CL); -%title(titles(ct)); -colormap(ax(2),'bone'); -xticks([]); -yticks([]); -end -plotWindow=[-0.5 2.5]; -plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); -for ct=1:8 -sel=category==ct; -activity=[]; -for cue=1:3 - for os=1:2 - activity = [activity PSTH3{cue,os}(sel,plotBins)]; - end -end -activity=activity./max(abs(activity),[],2); -[coeff,score,~,~,explained]=pca(activity'); -specs={'-','--'}; -for comp=1 - cond=0; - if explained(comp)>5 - subplot(4,8,ct) - - hold on - for cue=1:3 - for os=1:2 - cond=cond+1; - plot(binTimes(plotBins),-score((cond-1)*sum(plotBins)+1:cond*sum(plotBins),comp),specs{os},'color',colors{cue},'linewidth',1); - end - end - ylabel(sprintf('#%g (%g%%)',comp,round(explained(comp),1))); - yticks([]); - xticks([0 2.5]); - if comp==1 title(titles{ct}); end - -% plot([0 0],[min(score(:,comp)) max(score(:,comp))],':','color','k','linewidth',0.5); -% plot(plotWindow,[0 0],':','color','k','linewidth',0.5); - - end - xlim(cueWindow); end -end %% plot PC1 value from each region figure; @@ -2077,7 +2007,7 @@ plotWindow=[-0.5 2.5]; plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); for reg=1:length(regions) - sel=ismember(division,regions(reg)) & category==1; + sel=ismember(division,regions(reg)) & ismember(category,[2]); activity=[]; for cue=1:3 for os=1:2 @@ -2088,7 +2018,7 @@ [coeff,score,~,~,explained]=pca(activity'); specs={'-','--'}; comp=1; -if strcmp(regions(reg),'ACA') comp=2; end +if strcmp(regions(reg),'ACA') comp=1; end cond=0; if explained(comp)>5 subplot(4,9,reg) @@ -2118,7 +2048,7 @@ division=neuronGroup; for reg=1:4 - sel=ismember(division,reg) & category==1; + sel=ismember(division,reg) & category==2; activity=[]; for cue=1:3 for os=1:2 @@ -2155,228 +2085,7 @@ end -%% PCA value analysis for each region -figure; -numstraps=5000; -numn=10; -neuronRegionOlf=neuronRegionAdj; -neuronRegionOlf(ismember(neuronRegionAdj,{'EPd','PIR'}))={'OLF'}; -division=neuronRegionOlf; - -regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; -reggrps=[1:3]; - -plotWindow=[0 2.5]; -plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); -plotBinTimes=binTimes(plotBins); -meanWindow=[1 2.5]; -meanBins=plotBinTimes>=meanWindow(1) & plotBinTimes<=meanWindow(2); -baseWindow=[0 0.5]; -baseBins=plotBinTimes>=baseWindow(1) & plotBinTimes<=baseWindow(2); - -cats={1,2:8}; -valScoreStrap=cell(2); -for c=1:2 - valScoreStrap{c}=NaN(numstraps,length(regions)); - -for reg=1:length(regions) - sel=ismember(division,regions(reg)) & ismember(category,cats{c}); - activitySep=cell(2,1); - activity=[]; - for cue=1:3 - for os=1:2 - activity = [activity PSTH3{cue,os}(sel,plotBins)]; - activitySep{os} = [activitySep{os} PSTH3{cue,os}(sel,plotBins)]; - end - end - %activity=activity./max(abs(activity),[],2); - for os=1:2 - activitySep{os} = activitySep{os}./max(abs(activity),[],2); - end - - if strcmp(regions(reg),'PL') & c==1 - [coeff,projscore{1},~,~,explained]=pca(activitySep{1}'); - projscore{2}=activitySep{2}' * coeff; - - - specs={'-','--'}; - for os=1:2 - - subplot(2,2,os); - hold on; - cond=0; - for cue=1:3 - cond=cond+1; - plot3(projscore{os}((cond-1)*sum(plotBins)+1:cond*sum(plotBins),1),projscore{os}((cond-1)*sum(plotBins)+1:cond*sum(plotBins),2),projscore{os}((cond-1)*sum(plotBins)+1:cond*sum(plotBins),3),specs{os},'color',colors{cue},'linewidth',1); - end - if os==1 - ylabel('PC2'); - xlabel('PC1'); - zlabel('PC3'); - title(regions{reg}); - end - % yticks([]); - % xticks([]); - % zticks([]); - % axis([-4 15 -4 15 -4 15]); - view(160,60); - end - - - - end - - for strap=1:numstraps - entries=randsample(sum(sel),numn,'false'); - [coeff,projscore{1},~,~,explained]=pca(activitySep{1}(entries,:)'); - projscore{2}=activitySep{2}(entries,:)' * coeff; - - %dimension 1 is time - %dimension 2 is cue - %dimension 3 is component - vectors=reshape(projscore{2}(:,1:3),[sum(plotBins),3,3]); - - distPM=NaN(sum(plotBins),1); - distPF=NaN(sum(plotBins),1); - distFM=NaN(sum(plotBins),1); - distP=NaN(sum(plotBins),1); - distF=NaN(sum(plotBins),1); - distM=NaN(sum(plotBins),1); - baseline=mean([squeeze(mean(vectors(baseBins,1,:))),squeeze(mean(vectors(baseBins,2,:))),squeeze(mean(vectors(baseBins,3,:)))],2); - alldistances=pdist2(projscore{2}(:,1:3),projscore{2}(:,1:3),'euclidean'); - maxDistance=max(alldistances,[],'all'); - for bin=1:sum(plotBins) -% distPM(bin,1)=norm(squeeze(vectors(bin,1,:)-vectors(bin,3,:))); -% distPF(bin,1)=norm(squeeze(vectors(bin,1,:)-vectors(bin,2,:))); -% distFM(bin,1)=norm(squeeze(vectors(bin,2,:)-vectors(bin,3,:))); - distP(bin,1)=norm(squeeze(vectors(bin,1,:))-baseline)/maxDistance; - distF(bin,1)=norm(squeeze(vectors(bin,2,:))-baseline)/maxDistance; - distM(bin,1)=norm(squeeze(vectors(bin,3,:))-baseline)/maxDistance; - end - -% distPFM=distPF+distFM; -% -% discScore=distPM./maxDistance .* distPF./maxDistance .* distFM./maxDistance; -% axisScore=distPM ./ distPFM; -% valueScore=discScore.*axisScore; -% %valueScore=axisScore; - valScoreStrap{c,1}(strap,reg)=mean(distP(meanBins)); - valScoreStrap{c,2}(strap,reg)=mean(distF(meanBins)); - valScoreStrap{c,3}(strap,reg)=mean(distM(meanBins)); - end -end -end - - -for c=1:2 - subplot(4,2,4+c) -hold on -for cue=1:3 -errorbar([1:length(regions)]-0.2+0.2*(cue-1),mean(valScoreStrap{c,cue},1),std(valScoreStrap{c,cue},1),'o','color',colors{cue}); -end -title('Distance from baseline'); -ylim([0 1]); -xlim([0 length(regions)+1]); -xticks(1:length(regions)); -yticks([0 1]); - -end -bootpc=cell(3,2); -for c=1:2 -for cue=1:3 -bootpc{cue,c}=NaN(length(regions)); -for con1=1:length(regions) - for con2=1:length(regions) - bootpc{cue,c}(con1,con2)=(sum(sum(valScoreStrap{c,cue}(:,con1)<=valScoreStrap{c,cue}(:,con2)'))+1)/(numstraps^2+1); - end -end -end -end - -rankp=cell(1,2); -for c=1:2 -rankp{c}=NaN(2,length(regions)); -for reg=1:length(regions) - rankp{c}(1,reg)=(sum(sum(valScoreStrap{c,1}(:,reg)<=valScoreStrap{c,2}(:,reg)'))+1)/(numstraps^2+1); - rankp{c}(2,reg)=(sum(sum(valScoreStrap{c,2}(:,reg)<=valScoreStrap{c,3}(:,reg)'))+1)/(numstraps^2+1); -end -end - -%for region groups -numn=60; -valScoreStrap=cell(2); -for c=1:2 - valScoreStrap{c}=NaN(numstraps,length(reggrps)); - -for reg=1:length(reggrps) - sel=ismember(neuronGroup,reggrps(reg)) & ismember(category,cats{c}); - activitySep=cell(2,1); - activity=[]; - for cue=1:3 - for os=1:2 - activity = [activity PSTH3{cue,os}(sel,plotBins)]; - activitySep{os} = [activitySep{os} PSTH3{cue,os}(sel,plotBins)]; - end - end - %activity=activity./max(abs(activity),[],2); - for os=1:2 - activitySep{os} = activitySep{os}./max(abs(activity),[],2); - end - - for strap=1:numstraps - entries=randsample(sum(sel),numn,'false'); - [coeff,projscore{1},~,~,explained]=pca(activitySep{1}(entries,:)'); - projscore{2}=activitySep{2}(entries,:)' * coeff; - - %dimension 1 is time - %dimension 2 is cue - %dimension 3 is component - vectors=reshape(projscore{2}(:,1:3),[sum(plotBins),3,3]); - - distPM=NaN(sum(plotBins),1); - distPF=NaN(sum(plotBins),1); - distFM=NaN(sum(plotBins),1); - distP=NaN(sum(plotBins),1); - distF=NaN(sum(plotBins),1); - distM=NaN(sum(plotBins),1); - baseline=mean([squeeze(mean(vectors(baseBins,1,:))),squeeze(mean(vectors(baseBins,2,:))),squeeze(mean(vectors(baseBins,3,:)))],2); - alldistances=pdist2(projscore{2}(:,1:3),projscore{2}(:,1:3),'euclidean'); - maxDistance=max(alldistances,[],'all'); - for bin=1:sum(plotBins) -% distPM(bin,1)=norm(squeeze(vectors(bin,1,:)-vectors(bin,3,:))); -% distPF(bin,1)=norm(squeeze(vectors(bin,1,:)-vectors(bin,2,:))); -% distFM(bin,1)=norm(squeeze(vectors(bin,2,:)-vectors(bin,3,:))); - distP(bin,1)=norm(squeeze(vectors(bin,1,:))-baseline)/maxDistance; - distF(bin,1)=norm(squeeze(vectors(bin,2,:))-baseline)/maxDistance; - distM(bin,1)=norm(squeeze(vectors(bin,3,:))-baseline)/maxDistance; - end - -% distPFM=distPF+distFM; -% -% discScore=distPM./maxDistance .* distPF./maxDistance .* distFM./maxDistance; -% axisScore=distPM ./ distPFM; -% valueScore=discScore.*axisScore; -% %valueScore=axisScore; - valScoreStrap{c,1}(strap,reg)=mean(distP(meanBins)); - valScoreStrap{c,2}(strap,reg)=mean(distF(meanBins)); - valScoreStrap{c,3}(strap,reg)=mean(distM(meanBins)); - end -end -end - - -for c=1:2 - subplot(4,3,9+c) -hold on -for cue=1:3 -errorbar([1:length(reggrps)]-0.2+0.2*(cue-1),mean(valScoreStrap{c,cue},1),std(valScoreStrap{c,cue},1),'o','color',colors{cue}); -end -ylim([0 1]); -yticks([0 1]); -xlim([0 length(reggrps)+1]); -xticks(1:length(reggrps)); -end %% GLM spatial plotting, 3D figure; @@ -2384,7 +2093,7 @@ %all neuronXYZmm=neuronXYZ/1000; -cats={1,2:6,7}; +cats={1,2,3}; for ct=1:3 subplot(1,3,ct); hold on; @@ -2406,21 +2115,22 @@ end + %% coding dimension for CS+/CS50/CS- figure; %categories for subselection sels={}; -sels{1}=ismember(category,8); -sels{2}=ismember(category,2:6); +sels{1}=ismember(category,3); +sels{2}=ismember(category,2); sels{3}=category==1; colors{1,1}=[0.1 0.6 0.2]; colors{2,1}=[0.4 0.1 0.4]; colors{3,1}=[0.3 0.3 0.3]; -vcolors={[0.05 0.25 0.45],[0.7 0 1],[0 0.7 1]}; +vcolors={[0.6 0.6 0.6],[0.7 0 1],[0 0.7 1]}; numstraps=5000; amin=-0.1; -amax=1.2; +amax=1.5; cueWindow=[0 2.5]; cueBins=binTimes>=cueWindow(1) & binTimes<=cueWindow(2); @@ -2528,7 +2238,7 @@ subplot(3,4,2+(cl-1)*4); hold on; for cue=1:3 - plot(oproj{cue,1}(plotBins),oproj{cue,3}(plotBins),':','linewidth',1,'color',colors{cue}); +% plot(oproj{cue,1}(plotBins),oproj{cue,3}(plotBins),':','linewidth',1,'color',colors{cue}); plot(oproj{cue,1}(cueBins),oproj{cue,3}(cueBins),'linewidth',1,'color',colors{cue}); end @@ -2594,7 +2304,7 @@ ylabel('dist from cs-'); xlabel('frequency'); -legend('odor','trial type','value'); +legend('odor','meaning','value'); bootp=(sum(sum(baseCSMdist(:,1)<=baseCSMdist(:,3)'))+1)/(numstraps^2+1); text(0.1,0.9,sprintf('p = %g',round(bootp,2,'significant'))); bootp=(sum(sum(baseCSMdist(:,2)<=baseCSMdist(:,3)'))+1)/(numstraps^2+1); @@ -2609,15 +2319,15 @@ ylabel('CS+, CS50 angle'); xlabel('frequency'); -legend('odor','trial type','value'); +legend('odor','meaning','value'); bootp=(sum(sum(cspfangle(:,1)<=cspfangle(:,3)'))+1)/(numstraps^2+1); text(0.05,25,sprintf('p = %g',round(bootp,2,'significant'))); bootp=(sum(sum(cspfangle(:,2)<=cspfangle(:,3)'))+1)/(numstraps^2+1); text(0.15,25,sprintf('p = %g',round(bootp,2,'significant'))); bootp=(sum(sum(cspfangle(:,1)>=cspfangle(:,2)'))+1)/(numstraps^2+1); text(0.05,95,sprintf('p = %g',round(bootp,2,'significant'))); -ylim([0 90]); -yticks([0 45 90]); +ylim([0 45]); +yticks([0:15:45]); subplot(2,4,4); @@ -2628,7 +2338,7 @@ ylabel('same/other correlation'); xlabel('frequency'); -legend('odor','trial type','value'); +legend('odor','meaning','value'); bootp=(sum(sum(projCorrStrp(:,1)>=projCorrStrp(:,3)'))+1)/(numstraps^2+1); text(0.05,.97,sprintf('p = %g',round(bootp,3,'significant'))); bootp=(sum(sum(projCorrStrp(:,2)>=projCorrStrp(:,3)'))+1)/(numstraps^2+1); @@ -2638,321 +2348,232 @@ ylim([0.9 1]); yticks([0.9 0.95 1]); +%% perform value analysis with trial history +runglm=false; -%% perform regression for shuffling trial values for value cells -alreadyRun=true; - -if ~alreadyRun -tic - -shuffles=1000; %how many times to shuffle trial values within cue -components=10; folds=4; -NS=0; -varExpValueTrlSh=NaN(totalNeurons,1000); +varExpValueHist=NaN(totalNeurons,1); +%kernels=NaN(totalNeurons,31); binsPerCue=diff(windows{1})/binSize; -[~,bestValue]=max(varExpValueSh(:,3:4),[],2); -NN=0; -for session=1:sum(includedSessions) - - NS=NS+1; - - A = Xall{NS}; - rA = A * bR(:,1:components); - trains = getfolds(rA,folds,binsPerTrial); - - - %get original cue order for shuffle - firstBins={}; - cueOrder=[1 3 5 2 4 6]; - for cue=1:6 - firstBins{cue}=find(A(:,1+(cue-1)*binsPerCue)); - firstBins{cue}=[firstBins{cue} ones(length(firstBins{cue}),1)*cueOrder(cue)]; - end - alltrials=cat(1,firstBins{:}); - [sortedtrials,ordered]=sort(alltrials(:,1)); - originalOrder=alltrials(ordered,2); - - %make new shuffled predictor matrices - A3Sh={}; - for perm=1:size(varExpValueTrlSh,2) - - newOrder=originalOrder; - newValue=NaN(length(alltrials),1); - for cue=1:6 - originalValues=trialValues{session}(originalOrder==cue); - newValue(newOrder==cue)=originalValues(randsample(sum(originalOrder==cue),sum(newOrder==cue),'false')); - end +opts.alpha=0.5; +glmopts=glmnetSet(opts); - A3Sh{perm}=Xall3{session}; - for bin=1:binsPerCue - A3Sh{perm}(sortedtrials(:,1)-1+bin,bin)=newValue; - end - end - - for neuron=1:size(Yall{NS},2) - NN=NN+1; - if category(NN,1)==1 & bestValue(NN,1)==1 - y=Yall{NS}(:,neuron); - - - thisLambda=lambdaVal(NN,1); - for perm=1:size(varExpValueTrlSh,2) - - rA3 = A3Sh{perm} * bR3(:,1:components); - - % if ismember(perm,topModels) - % display(perms6all(perm,:));display([combinedValues(1:10) newValue(1:10)]); - % end - - - pred=NaN(size(y,1),1); - for fold=1:folds - train=trains{fold}; - test=train==0; - fitK=lassoglm(rA3(train,:),y(train),'normal','alpha',0.5,'lambda',thisLambda); - kernel=bR3(:,1:components)*fitK; - pred(test)=A3Sh{perm}(test,:)*kernel; - end - varExpValueTrlSh(NN,perm) = 1- var(y-pred)/var(y); - end - end - +if runglm +NN=0; +sessInd=find(includedSessions); +tic +for NS=1:sum(includedSessions) + A = Xall3{NS}(:,[1:binsPerCue size(Xall3{NS},2)-5:size(Xall3{NS},2)]); + trains = getfolds(A,folds,binsPerTrial); + NN=NNstart; + for neuron=1:size(Yall{NS},2) + NN=NN+1; + y=Yall{NS}(:,neuron); + %cross-validated variance explained + pred=NaN(size(y,1),1); + for fold=1:folds + train=trains{fold}; + test=train==0; + %kernel=lassoglm(A(train,:),y(train),'normal','alpha',0.5,'lambda',thisLambda); + fit = glmnet(A(train,:),y(train),[],glmopts); + prediction = glmnetPredict(fit,A(test,:)); + pred(test)=prediction(:,end); end - fprintf('Session #%d \n',NS); - -end - -save('newGLMFile.mat','varExp','varExpValueSh','varExpValueTrlSh','kernels','lambdaVal','predF','-v7.3'); - - -toc -end -%% history coding figure -figure; -%categories for subselection - -[~,bestValue]=max(varExpValueSh(:,3:4),[],2); -%how many shuffles was the history model better than? -overShuffle=sum(varExpValueSh(:,3)>varExpValueTrlSh,2)/size(varExpValueTrlSh,2); -%remove neurons that did poorly (better than fewer than 65% of shuffles) -bestValue(category==1&bestValue==1&overShuffle<0.65)=2; - -regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; -reggrps=[1:3]; - -% regions={'ACA','FRP','PL','ILA','ORB','CP','ACB'}; -% reggrps=[2 4]; - -allReg=unique(neuronRegionOlf); -regInd=NaN(length(regions),1); -for reg=1:length(regions) - regInd(reg,1)=find(ismember(allReg,regions(reg))); -end - - -for ct=1 - %region - barData=zeros(2,5); - subplot(2,2,1); - hold on; - sel=category==1 & bestValue==1; - othervalue=ismember(category,1)&~sel; - othersel=ismember(category,2:8); - regionnumber=zeros(totalNeurons,1); - for reg=1:length(allReg) - regsel=ismember(neuronRegionOlf,allReg(reg)); - regionnumber(regsel)=reg; - barData(1,reg)=sum(sel®sel)/sum(regsel); - barData(2,reg)=sum(othervalue®sel)/sum(regsel); - barData(3,reg)=sum(othersel®sel)/sum(regsel); - end - tbl = table(categorical(neuronSession(regionnumber>0)), categorical(regionnumber(regionnumber>0)), sel(regionnumber>0), 'VariableNames', {'session', 'region','val'}); - lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial'); - [ypred,ypredCI]=predict(lm,'conditional',false); - - barData=barData(:,regInd); - b=bar(barData','stacked'); - b(2).FaceColor=[0 0.7 1]; - b(3).FaceColor=[0.85 0.85 0.85]; - b(1).FaceColor=[0.4 0.9 1]; - - upperCI=[]; - lowerCI=[]; - estmean=[]; - for reg=1:length(regions) - regsel=find(ismember(neuronRegionOlf,regions(reg)),1); - estmean(1,reg)=ypred(regsel); - upperCI(1,reg)=ypredCI(regsel,2); - lowerCI(1,reg)=ypredCI(regsel,1); - end - - errorbar(1:length(regions),estmean(1,:),estmean(1,:)-lowerCI(1,:),upperCI(1,:)-estmean(1,:),'o','linewidth',0.75,'color','k'); - - xlim([0.5 length(regions)+0.5]); - %legend('Cues','Cues,Licks','Licks','Licks,Reward','Reward','Cues,Reward','Cues,Licks,Reward'); - xticks(1:length(regions)); - xticklabels(regions); - xtickangle(45); - ylim([0 0.55]); - - - %make heatmap - target=[0.4 0.9 1]; - customMap=NaN(100,3); - for i=1:length(customMap) - customMap(i,:)=target*i/(length(customMap)); - end - customMap=cat(1,[1 1 1],customMap); - - %by region - greater=zeros(length(regions)); - frac=zeros(length(regions)); - for r1=1:length(regions) - for r2=1:length(regions) - frac(r1,r2)=barData(1,r1)-barData(1,r2); - greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); - end - end - s=subplot(2,2,2); - imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); - - colormap(s,customMap); - cb=colorbar; - set(cb,'ticks',[0:0.05:0.3]); - xticks(1:length(regions)); - xticklabels(regions); - xtickangle(45); - yticks(1:length(regions)); - yticklabels(regions); - %title(titles{ct}); - if ct==1 - ylabel('increase in this region...'); - xlabel('over this region'); + varExpValueHist(NN,1) = 1- var(y-pred)/var(y); + kernel = fit.beta(:,end); + kernels(NN,:)=kernel; end + NNstart=NN; + fprintf('Session #%d \n',NS); +end +toc +save('valueGLMFile.mat','varExp','varExpValueSh','varExpValueNew','varExpValueHist','kernels','lambdaVal','predF','-v7.3'); + +end + + +%% perform regression for shuffling trial values for value cells +alreadyRun=true; + +if ~alreadyRun +tic + +shuffles=1000; %how many times to shuffle trial values within cue +folds=4; +NS=0; +varExpValueTrlSh=NaN(totalNeurons,shuffles); +binsPerCue=diff(windows{1})/binSize; +NN=0; + +for session=1:sum(includedSessions) - %region group - barData=zeros(2,3); - subplot(2,2,3); - hold on; - regionnumber=zeros(totalNeurons,1); - for reg=1:4 - regsel=ismember(neuronGroup,reg); - regionnumber(regsel)=reg; - barData(1,reg)=sum(sel®sel)/sum(regsel); - barData(2,reg)=sum(othervalue®sel)/sum(regsel); - barData(3,reg)=sum(othersel®sel)/sum(regsel); - end - tbl = table(categorical(neuronSession(regionnumber>0)), categorical(regionnumber(regionnumber>0)), sel(regionnumber>0), 'VariableNames', {'session', 'region','val'}); - lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial'); - [ypred,ypredCI]=predict(lm,'conditional',false); - - barData=barData(:,reggrps); - b=bar(barData','stacked'); - b(2).FaceColor=[0 0.7 1]; - b(3).FaceColor=[0.85 0.85 0.85]; - b(1).FaceColor=[0.4 0.9 1]; + NS=NS+1; - upperCI=[]; - lowerCI=[]; - estmean=[]; - for reg=1:length(reggrps) - regsel=find(ismember(neuronGroup,reggrps(reg)),1); - estmean(1,reg)=ypred(regsel); - upperCI(1,reg)=ypredCI(regsel,2); - lowerCI(1,reg)=ypredCI(regsel,1); + A = Xall{NS}; + trains = getfolds(A,folds,binsPerTrial); + + %get original cue order for shuffle + firstBins={}; + cueOrder=[1 3 5 2 4 6]; + for cue=1:6 + firstBins{cue}=find(A(:,1+(cue-1)*binsPerCue)); + firstBins{cue}=[firstBins{cue} ones(length(firstBins{cue}),1)*cueOrder(cue)]; end + alltrials=cat(1,firstBins{:}); + [sortedtrials,ordered]=sort(alltrials(:,1)); + originalOrder=alltrials(ordered,2); - errorbar(1:length(reggrps),estmean(1,:),estmean(1,:)-lowerCI(1,:),upperCI(1,:)-estmean(1,:),'o','linewidth',0.75,'color','k'); - - xlim([0.5 7+0.5]); - %legend('Cues','Cues,Licks','Licks','Licks,Reward','Reward','Cues,Reward','Cues,Licks,Reward'); - xticks(1:length(reggrps)); - xticklabels(regionGroupNames(reggrps)); - xtickangle(45); - ylim([0 0.55]); - - %by region - greater=zeros(length(reggrps)); - frac=zeros(length(reggrps)); - for r1=1:length(reggrps) - for r2=1:length(reggrps) - frac(r1,r2)=barData(1,r1)-barData(1,r2); - greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); - end + %make new shuffled predictor matrices + A3Sh={}; + for perm=1:size(varExpValueTrlSh,2) + + newOrder=originalOrder; + newValue=NaN(length(alltrials),1); + for cue=1:6 + originalValues=trialValues{session}(originalOrder==cue); + newValue(newOrder==cue)=originalValues(randsample(sum(originalOrder==cue),sum(newOrder==cue),'false')); + end + + A3Sh{perm}=Xall3{session}(:,[1:binsPerCue size(Xall3{NS},2)-5:size(Xall3{NS},2)]); + for bin=1:binsPerCue + A3Sh{perm}(sortedtrials(:,1)-1+bin,bin)=newValue; + end end - s=subplot(2,2,4); - imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); - colormap(s,customMap); - cb=colorbar; - set(cb,'ticks',[0:0.05:0.3]); - xticks(1:length(reggrps)); - xticklabels(regionGroupNames(reggrps)); - xtickangle(45); - yticks(1:length(reggrps)); - yticklabels(regionGroupNames(reggrps)); - %title(titles{ct}); - if ct==1 - ylabel('increase in this region...'); - xlabel('over this region'); - end + for neuron=1:size(Yall{NS},2) + NN=NN+1; + if ismember(category(NN,1),[1 2]) + + y=Yall{NS}(:,neuron); + + parfor perm=1:size(varExpValueTrlSh,2) + A=A3Sh{perm}; + %cross-validated variance explained + pred=NaN(size(y,1),1); + for fold=1:folds + train=trains{fold}; + test=train==0; + %kernel=lassoglm(A(train,:),y(train),'normal','alpha',0.5,'lambda',thisLambda); + fit = glmnet(A(train,:),y(train),[],glmopts); + prediction = glmnetPredict(fit,A(test,:)); + pred(test)=prediction(:,end); + end + varExpValueTrlSh(NN,perm) = 1- var(y-pred)/var(y); + + + + end + + end + end + fprintf('Session #%d \n',NS); end +save('valueGLMFile.mat','varExp','varExpValueSh','varExpValueNew','varExpValueHist','kernels','varExpValueTrlSh','predF','-v7.3'); +toc +end +%% history coding figure +figure; +%categories for subselection +sel = ismember(category,[1:8]); %which neurons are allowed +%usedModels=varExpValueNew(:,91:end); +models2compare = [varExpValueHist varExpValueNew(:,91:end)]; +%models2compare = [usedModels(:,valRank(1:17))]; +[~,bestModelH]=max(models2compare,[],2); +overShuffle=sum(varExpValueHist>varExpValueTrlSh,2)/size(varExpValueTrlSh,2); +history=bestModelH==1 & overShuffle>0.95 & sel; +neuronRegionOlf=neuronRegionAdj; +neuronRegionOlf(ismember(neuronRegionAdj,{'EPd','PIR'}))={'OLF'}; +allReg=unique(neuronRegionOlf); - -%% history coding proportions out of value cells regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; -reggrps=[1:3]; +reggrps=[2:4]; +neuronGroupAdj = neuronGroup + 1; +neuronGroupAdj(neuronGroupAdj==5)=1; -%regions={'ACA','FRP','PL','ILA','ORB','CP','ACB'}; -%reggrps=[2 4]; +% regions={'ACA','FRP','PL','ILA','ORB','CP','ACB'}; +% reggrps=[2 4]; allReg=unique(neuronRegionOlf); regInd=NaN(length(regions),1); +order=[9 10 5 2 3 1 7 12 8 11 6 13 4]; +allReg=allReg(order); for reg=1:length(regions) regInd(reg,1)=find(ismember(allReg,regions(reg))); end +pVals=NaN(36,1); for ct=1 %region barData=zeros(2,5); subplot(2,2,1); hold on; - sel=category==1 & bestValue==1; - othervalue=ismember(category,1)&~sel; - othersel=ismember(category,2:8); + sel=history; + %othervalue=ismember(category,1)&~sel; + othersel=ismember(category,1:8)&~sel; regionnumber=zeros(totalNeurons,1); for reg=1:length(allReg) - regsel=ismember(neuronRegionOlf,allReg(reg))&category==1; + regsel=ismember(neuronRegionOlf,allReg(reg)); regionnumber(regsel)=reg; barData(1,reg)=sum(sel®sel)/sum(regsel); - barData(2,reg)=sum(othervalue®sel)/sum(regsel); - barData(3,reg)=sum(othersel®sel)/sum(regsel); + %barData(2,reg)=sum(othervalue®sel)/sum(regsel); + barData(2,reg)=sum(othersel®sel)/sum(regsel); end tbl = table(categorical(neuronSession(regionnumber>0)), categorical(regionnumber(regionnumber>0)), sel(regionnumber>0), 'VariableNames', {'session', 'region','val'}); - lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial'); + lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial','dummyvarcoding','reference'); [ypred,ypredCI]=predict(lm,'conditional',false); - - barData=barData(:,regInd); - b=bar(barData','stacked'); - b(2).FaceColor=[0 0.7 1]; - b(3).FaceColor=[0.85 0.85 0.85]; - b(1).FaceColor=[0.4 0.9 1]; - + + Hall=[]; +% Htemp = [1 0 0 0 0 0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0 0 0 0 0 0]; +% for reg1 = 1:length(regions) +% H=Htemp; +% H(2,reg1+4)=1; +% Hall=cat(1,Hall,H); +% end + Htemp = [0 0 0 0 0 0 0 0 0 0 0 0 0]; + for reg1 = 1:length(regions)-1 + for reg2 = reg1+1:length(regions) + H=Htemp; + H(reg1+4)=1; + H(reg2+4)=-1; + Hall=cat(1,Hall,H); + end + end + greater=zeros(length(regions)); + for reg1 = 1:length(regions) + for reg2 = 1:length(regions) + H=Htemp; + H(reg1+4)=1; + H(reg2+4)=-1; + greater(reg1,reg2)=coefTest(lm,H); + end + end + + %Hall=unique(Hall,'rows'); + + for p=1:length(pVals) + pVals(p,ct) = coefTest(lm,Hall(p,:)); + end + upperCI=[]; lowerCI=[]; estmean=[]; for reg=1:length(regions) - regsel=find(ismember(neuronRegionOlf(category==1),regions(reg)),1); + regsel=find(ismember(neuronRegionOlf,regions(reg)),1); estmean(1,reg)=ypred(regsel); upperCI(1,reg)=ypredCI(regsel,2); lowerCI(1,reg)=ypredCI(regsel,1); - end + end + + barData=barData(:,regInd); + b=bar(barData','stacked'); + b(2).FaceColor=[0.85 0.85 0.85]; %meaning + b(1).FaceColor=[0.4 0.9 1]; %value + errorbar(1:length(regions),estmean(1,:),estmean(1,:)-lowerCI(1,:),upperCI(1,:)-estmean(1,:),'o','linewidth',0.75,'color','k'); @@ -2961,7 +2582,7 @@ xticks(1:length(regions)); xticklabels(regions); xtickangle(45); - ylim([0 1]); + ylim([0 0.55]); %make heatmap @@ -2973,16 +2594,17 @@ customMap=cat(1,[1 1 1],customMap); %by region - greater=zeros(length(regions)); + %greater=zeros(length(regions)); frac=zeros(length(regions)); for r1=1:length(regions) for r2=1:length(regions) frac(r1,r2)=barData(1,r1)-barData(1,r2); - greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); + %greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end s=subplot(2,2,2); - imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); + pcutoff=0.05/36; + imagesc(frac .* (greater0), [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); cb=colorbar; @@ -3004,11 +2626,11 @@ hold on; regionnumber=zeros(totalNeurons,1); for reg=1:4 - regsel=ismember(neuronGroup,reg)&category==1; + regsel=ismember(neuronGroupAdj,reg); regionnumber(regsel)=reg; barData(1,reg)=sum(sel®sel)/sum(regsel); - barData(2,reg)=sum(othervalue®sel)/sum(regsel); - barData(3,reg)=sum(othersel®sel)/sum(regsel); + %barData(2,reg)=sum(othervalue®sel)/sum(regsel); + barData(2,reg)=sum(othersel®sel)/sum(regsel); end tbl = table(categorical(neuronSession(regionnumber>0)), categorical(regionnumber(regionnumber>0)), sel(regionnumber>0), 'VariableNames', {'session', 'region','val'}); lm = fitglme(tbl,'val~1+region+(1|session)','distribution','binomial'); @@ -3016,15 +2638,15 @@ barData=barData(:,reggrps); b=bar(barData','stacked'); - b(2).FaceColor=[0 0.7 1]; - b(3).FaceColor=[0.85 0.85 0.85]; + %b(2).FaceColor=[0 0.7 1]; + b(2).FaceColor=[0.85 0.85 0.85]; b(1).FaceColor=[0.4 0.9 1]; upperCI=[]; lowerCI=[]; estmean=[]; for reg=1:length(reggrps) - regsel=find(ismember(neuronGroup(category==1),reggrps(reg)),1); + regsel=find(ismember(neuronGroupAdj,reggrps(reg)),1); estmean(1,reg)=ypred(regsel); upperCI(1,reg)=ypredCI(regsel,2); lowerCI(1,reg)=ypredCI(regsel,1); @@ -3035,30 +2657,43 @@ xlim([0.5 7+0.5]); %legend('Cues','Cues,Licks','Licks','Licks,Reward','Reward','Cues,Reward','Cues,Licks,Reward'); xticks(1:length(reggrps)); - xticklabels(regionGroupNames(reggrps)); + xticklabels(regionGroupNames(reggrps-1)); xtickangle(45); - ylim([0 1]); + ylim([0 0.55]); + - %by region + Htemp = [0 0 0 0]; greater=zeros(length(reggrps)); + for reg1 = 1:length(reggrps) + for reg2 = 1:length(reggrps) + H=Htemp; + H(reg1+1)=1; + H(reg2+1)=-1; + greater(reg1,reg2)=coefTest(lm,H); + end + end + + %by region + %greater=zeros(length(reggrps)); frac=zeros(length(reggrps)); for r1=1:length(reggrps) for r2=1:length(reggrps) frac(r1,r2)=barData(1,r1)-barData(1,r2); - greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); + %greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end s=subplot(2,2,4); - imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); + pcutoff=0.05/3; + imagesc(frac .* (greater0), [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); cb=colorbar; set(cb,'ticks',[0:0.05:0.3]); xticks(1:length(reggrps)); - xticklabels(regionGroupNames(reggrps)); + xticklabels(regionGroupNames(reggrps-1)); xtickangle(45); yticks(1:length(reggrps)); - yticklabels(regionGroupNames(reggrps)); + yticklabels(regionGroupNames(reggrps-1)); %title(titles{ct}); if ct==1 ylabel('increase in this region...'); @@ -3066,17 +2701,18 @@ end end -%% value models schematic -examplePerms=[1 1]; -exampleNeurons=[5 20 50 72 111 141]; -%5 12 20 36 37 50 72 111 141 161 + +%% history value models schematic +examplePerms=[1 1 1]; +exampleNeurons=[31 37 63]; +%31 37 63 plotWindow=[0 2.5]; plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); binTimesGLMTrl = analysisWindow(1)+binSize/2:binSize:analysisWindow(2)-binSize/2; plotBinsGLM=binTimesGLMTrl>=plotWindow(1) & binTimesGLMTrl<=plotWindow(2); - +cueValues=[1 1 0.5 0.5 0 0]; valcolors{1,1}=[0.2 0.2 0.2]; valcolors{2,1}=[0.5 0.5 0.5]; valcolors{3,1}=[0.35 0.05 0.35]; @@ -3088,17 +2724,11 @@ valcolors{9,1}=[0.1 0.6 0.2]; valcolors{10,1}=[0.2 0.8 0.3]; -%smoothing filter for licks -smoothbinsLick=25; %number of previous bins used to smooth -halfnormal=makedist('HalfNormal','mu',0,'sigma',20); -filterweightsLick=pdf(halfnormal,0:smoothbinsLick); - -components=10; folds=4; neuronNumber=[1:length(neuronSession)]'; %get neuron number from order in value category -entries=find(category==1 & bestValue==1); +entries=find(history); for n=1:length(exampleNeurons) nn=entries(exampleNeurons(n)); ta=[]; @@ -3128,114 +2758,50 @@ session=neuronSession(nn); sessionNeuron=find(neuronNumber(neuronSession==session)==nn); sessionFolder=fullfile(direc,sessionSubject{session},sessionDate{session}); - lickTimes = readNPY(fullfile(sessionFolder,'licks.times.npy')); - reward = readNPY(fullfile(sessionFolder,'rewards.times.npy')); cueTimes = readNPY(fullfile(sessionFolder,'cues.times.npy')); rewprob = readNPY(fullfile(sessionFolder,'cues.rewardProbability.npy')); odorset = readNPY(fullfile(sessionFolder,'cues.odorSet.npy')); + cueOrder=zeros(length(cueTimes),1); + cueOrder(rewprob==1&odorset==1)=1; + cueOrder(rewprob==1&odorset==2)=2; + cueOrder(rewprob==0.5&odorset==1)=3; + cueOrder(rewprob==0.5&odorset==2)=4; + cueOrder(rewprob==0&odorset==1)=5; + cueOrder(rewprob==0&odorset==2)=6; csp1 = cueTimes(rewprob==1&odorset==1); csp2 = cueTimes(rewprob==1&odorset==2); csf1 = cueTimes(rewprob==0.5&odorset==1); csf2 = cueTimes(rewprob==0.5&odorset==2); csm1 = cueTimes(rewprob==0&odorset==1); - csm2 = cueTimes(rewprob==0&odorset==2); - - %first lick in bout - ibi=0.5; - boutStart = lickTimes(cat(1,1,1+find(diff(lickTimes)>ibi))); - - binEdgesSession=cueTimes(1)-5:binSizeSession:cueTimes(end)+10; - binTimesSession=cueTimes(1)-5+binSizeSession/2:binSizeSession:cueTimes(end)+10-binSizeSession/2; - binnedLicksGLM=histcounts(lickTimes,binEdgesSession) / binSizeSession; - lickRateGLM=NaN(length(binnedLicksGLM),1); - for l=1:length(binnedLicksGLM) - lickRateGLM(l,1)=sum(binnedLicksGLM(l-min([l-1 smoothbinsLick]):l).*fliplr(filterweightsLick(1:min([l smoothbinsLick+1]))))/sum(filterweightsLick(1:min([l smoothbinsLick+1]))); - end - - binnedLicksGLM=NaN(length(cueTimes),length(binEdges)-1); - binTimesGLM=NaN(length(cueTimes),length(binEdges)-1); - for trial = 1:length(cueTimes) - binTimesRel = (binTimesSession - cueTimes(trial)); - for bin=1:length(binTimes) - binnedLicksGLM(trial,bin)=mean(lickRateGLM(binTimesRel>=binEdges(bin)&binTimesRelanalysisWindow(1) & binTimes=cueWindow(1) & binTimes<=cueWindow(2); @@ -3380,6 +2946,7 @@ end bootpp={}; +valSlope=NaN(numstraps,length(sels)); for cl=1:length(sels) sel=sels{cl}; @@ -3414,7 +2981,7 @@ end xticks([0 1 2]); yticks([0 1]); - plot([0 0],[-2 12],':','color','k','linewidth',0.5); + plot([0 0],[-3 12],':','color','k','linewidth',0.5); plot(plotWindow,[0 0],':','color','k','linewidth',0.5); for cue=1:6 @@ -3452,7 +3019,8 @@ position(strap,cue)=mean(sproj{cue,1}); end - + b=[(((.225:0.05:.475)-0.05)/0.45);ones(1,6)]'\position(strap,3:8)'; + valSlope(strap,cl)=b(1); end @@ -3463,6 +3031,7 @@ end end + subplot(3,5,10+cl); hold on; for cue=1:c @@ -3481,14 +3050,118 @@ end -%% get proportions of value and trial type and compare between regions for striatum +catcolors={[0.4 0.9 1],[0 0.7 1],[0.7 0 1],[0.6 0.6 0.6],[0.9 0.2 0]}; +subplot(3,5,1); +hold on +for s=1:length(sels) + errorbar(s,nanmean(valSlope(:,s)),nanstd(valSlope(:,s)),'o','color',catcolors{s},'linewidth',1); +end +bootsp=NaN(length(sels)); +for con1=1:length(sels) + for con2=1:length(sels) + bootsp(con1,con2)=(sum(sum(valSlope(:,con1)<=valSlope(:,con2)'))+1)/(numstraps^2+1); + end +end +xticks(1:s); +yticks(0:0.5:1); +ylim([0 1]); +xlim([0.5 s+0.5]); +ylabel('CS50 value slope'); +basep=(sum(valSlope<=0,1)+1)/(numstraps+1); + +subplot(3,4,2); +hold on +plot([-0.1 1],[-0.1 1],':','linewidth',0.75,'color','k'); +scatter(varExpValueHist(history),varExp(history,2),24,catcolors{1},'filled');alpha(0.5); +%histogram(varExpValueHist(history)-varExp(history,2),'facecolor',catcolors{1}); +xlabel('Var. exp. with history model'); +ylabel('Var. exp. with lick model'); +xlim([-0.1 1]); +ylim([-0.1 1]); +xticks(0:0.5:1); +yticks(0:0.5:1); + +%% 1000 value shuffles +runglm=true; +%generate all value shuffles + +folds=4; +NNstart=0; +%varExpValueNew=NaN(totalNeurons,length(valueShuffles)); +%kernels=NaN(totalNeurons,31); +binsPerCue=diff(windows{1})/binSize; + +opts.alpha=0.5; +glmopts=glmnetSet(opts); + +if runglm +varExpValueShuffles=NaN(totalNeurons,1000); +sessInd=find(includedSessions); +valOrder=[1 1 0.5 0.5 0 0]; +tic +for NS=1:sum(includedSessions) + session=sessInd(incses); + sessionFolder=fullfile(direc,sessionSubject{session},sessionDate{session}); + cueTimes = readNPY(fullfile(sessionFolder,'cues.times.npy')); + rewprob = readNPY(fullfile(sessionFolder,'cues.rewardProbability.npy')); + odorset = readNPY(fullfile(sessionFolder,'cues.odorSet.npy')); + cueOrder=zeros(length(cueTimes),1); + cueOrder(rewprob==1&odorset==1)=1; + cueOrder(rewprob==1&odorset==2)=2; + cueOrder(rewprob==0.5&odorset==1)=3; + cueOrder(rewprob==0.5&odorset==2)=4; + cueOrder(rewprob==0&odorset==1)=5; + cueOrder(rewprob==0&odorset==2)=6; + + A = Xall{NS}(:,[1:binsPerCue size(Xall{NS},2)-5:size(Xall{NS},2)]); + trains = getfolds(A,folds,binsPerTrial); + for sh=1:1000 + cueOrderSh=cueOrder(randperm(length(cueOrder)),1); + for c=1:length(cueTimes) + cueValue=valOrder(cueOrderSh(c)); + dm=cueValue*eye(binsPerCue); + A(c*binsPerTrial-binsPerCue+1:c*binsPerTrial,1:binsPerCue)=dm; + end + NN=NNstart; + for neuron=1:size(Yall{NS},2) + NN=NN+1; + if category(NN)==1 + y=Yall{NS}(:,neuron); + %cross-validated variance explained + pred=NaN(size(y,1),1); + for fold=1:folds + train=trains{fold}; + test=train==0; + %kernel=lassoglm(A(train,:),y(train),'normal','alpha',0.5,'lambda',thisLambda); + fit = glmnet(A(train,:),y(train),[],glmopts); + prediction = glmnetPredict(fit,A(test,:)); + pred(test)=prediction(:,end); + end + varExpValueShuffles(NN,sh) = 1- var(y-pred)/var(y); + end + + end + end + fprintf('Session #%d \n',NS); + NNstart=NN; +end +toc +save('valueGLMFile.mat','varExp','varExpValueSh','varExpValueNew','varExpValueShuffles','varExpValueHist','kernels','varExpValueTrlSh','predF','-v7.3'); + +end + +%% how many neurons better than chance fitting of value? + +overShuffle=sum(varExpValueShuffles<=varExpValueNew(:,1),2)/size(varExpValueShuffles,2); +overShuffle=overShuffle(category==1); +%% get proportions of value and meaning and compare between regions for striatum regions={'ACA','FRP','PL','ILA','ORB','CP','ACB'}; %fit all together, allowing random effect of session to have more power catcolors={[0 0.7 1],[0.7 0 1],[0.6 0.6 0.6]}; -titles={'value','trial type (non-value)','non-trial type'}; -cats={1,2:6,7}; +titles={'value','value-like','untuned'}; +cats={1,2,3}; for ct=1:3 %region barData=zeros(2,5); @@ -3508,8 +3181,12 @@ [ypred,ypredCI]=predict(lm,'conditional',false); b=bar(barData','stacked'); - b(2).FaceColor=[0.85 0.85 0.85]; - b(1).FaceColor=catcolors{ct}; + %b(3).FaceColor=[0.4 0 0.5]; %odor + b(2).FaceColor=[0.85 0.85 0.85]; %meaning + b(1).FaceColor=catcolors{ct}; %value + %b(4).FaceColor=[0 0 0.2]; %odor + % b(5).FaceColor=[0.9 0.2 0]; + %b(3).FaceColor=[1 1 1]; upperCI=[]; lowerCI=[]; estmean=[]; diff --git a/Code/ephysValueDecoding.m b/Code/ephysValueDecoding.m new file mode 100644 index 0000000..9437708 --- /dev/null +++ b/Code/ephysValueDecoding.m @@ -0,0 +1,751 @@ +%% get raw spikes at different time bins for decoding +dwindow=[-0.5 2.5]; +dbinsize=0.25; +dbinEdges=dwindow(1):dbinsize:dwindow(2); +dbinTimes=dwindow(1)+dbinsize/2:dbinsize:dwindow(2)-dbinsize/2; + +preRun = true; + +if preRun + load('dactivity.mat'); +else + dactivity=cell(totalNeurons,6); + NN=0; + for incses=1:sum(includedSessions) + session=sessInd(incses); + %get events + sessionFolder=fullfile(direc,sessionSubject{session},sessionDate{session}); + cueTimes = readNPY(fullfile(sessionFolder,'cues.times.npy')); + rewprob = readNPY(fullfile(sessionFolder,'cues.rewardProbability.npy')); + odorset = readNPY(fullfile(sessionFolder,'cues.odorSet.npy')); + csp1 = cueTimes(rewprob==1&odorset==1); + csp2 = cueTimes(rewprob==1&odorset==2); + csf1 = cueTimes(rewprob==0.5&odorset==1); + csf2 = cueTimes(rewprob==0.5&odorset==2); + csm1 = cueTimes(rewprob==0&odorset==1); + csm2 = cueTimes(rewprob==0&odorset==2); + trialNo=1:length(cueTimes); + + folders=dir(sessionFolder); + probeFolderIDs = find(arrayfun(@(x) strcmp('p',folders(x).name(1)),1:length(folders))); + if isreal(probeFolderIDs) + for probe=1:length(probeFolderIDs) + probeFolder=fullfile(sessionFolder,folders(probeFolderIDs(probe)).name); + shankDir=dir(probeFolder); + entryName={}; + entryNameWithShank={}; + for entry=1:length(shankDir) + entryName{entry,1}=shankDir(entry).name(1:min([length(shankDir(entry).name) 5])); + entryNameWithShank{entry,1}=shankDir(entry).name(1:min([length(shankDir(entry).name) 6])); + end + shanks=sum(strcmp(entryName,'shank')); + shankList=[]; + for entry=1:length(shankDir) + if length(entryNameWithShank{entry})>5 + shankList=cat(1,shankList,str2num(entryNameWithShank{entry}(6))); + end + end + for shank=1:shanks + + shankFolder=fullfile(probeFolder,append('shank',num2str(shankList(shank)))); + clusterNames=readNPY(fullfile(shankFolder,'clusters.names.npy')); + incClusters=clusterNames(readNPY(fullfile(shankFolder,'clusters.passedQualityControl.npy'))); + channelLocations = textscan(fopen(fullfile(shankFolder,'channels.CCFregions.txt')),'%s'); + fclose('all'); + channelLocations = channelLocations{1}; + spikeTimes = readNPY(fullfile(shankFolder,'spikes.times.npy')); + spikeClusters = readNPY(fullfile(shankFolder,'spikes.clusterNames.npy')); + clusterChannel = readNPY(fullfile(shankFolder,'clusters.channels.npy')); + + %get clusters + for clus=1:length(incClusters) + + cluster = incClusters(clus); + regionName = channelLocations{clusterChannel(clusterNames==cluster)}; + charOnly = regexp(regionName,'[c-zA-Z]','match'); + noNums=cat(2,charOnly{:}); + if sum(strcmp(allAreas,noNums(1:min([3 length(noNums)]))))==1 + NN=NN+1; + clusterTimes = double(spikeTimes(spikeClusters==cluster))/30000; + + + binnedActivity=NaN(length(cueTimes),length(dbinEdges)-1); + numBins=length(dbinTimes); + for trial = 1:length(cueTimes) + clusterTimesRel = clusterTimes - cueTimes(trial); + binnedActivity(trial,:)=histcounts(clusterTimesRel,dbinEdges); + end + binnedActivity(isnan(binnedActivity))=0;binnedActivity(binnedActivity==inf)=0; %this only happens if bins go beyond end of session, so extremely rare + + %3 conditions + selections={}; + selections{1,1}=csp1; %cs+ + selections{2,1}=csp2; %cs+ + selections{3,1}=csf1; %cs50 + selections{4,1}=csf2; %cs50 + selections{5,1}=csm1; %cs- + selections{6,1}=csm2; %cs- + for condition=1:length(selections) + sel = ismember(round(cueTimes),round(selections{condition})); + dactivity{NN,condition}=binnedActivity(sel,:); + end + end + end + end + end + end + + + fprintf('Session %d \n',session); + end + + save('dactivity.mat','dactivity','-v7.3'); + +end + + +%% single unit decoding of cue identity +preRun=true; + +if preRun==false +dts=50; +folds=5; +trains=ones(dts*6,folds); +for t=1:folds + trains(t:5:end,t)=0; +end + +odorLabels=[ones(dts,1);2*ones(dts,1);3*ones(dts,1);4*ones(dts,1);5*ones(dts,1);6*ones(dts,1)]; +valueLabels=[ones(dts,1);ones(dts,1);0.5*ones(dts,1);0.5*ones(dts,1);zeros(dts,1);zeros(dts,1)]; +setLabels=[ones(dts,1);2*ones(dts,1);ones(dts,1);2*ones(dts,1);ones(dts,1);2*ones(dts,1)]; + +tps=size(dactivity{1,1},2); +odorPerf=NaN(totalNeurons,tps); +valuePerf=NaN(totalNeurons,tps); +setPerf=NaN(totalNeurons,tps); +tic +for NN=1:totalNeurons + + activity=[]; + for cue=1:6 + trls=randperm(length(dactivity{NN,cue}),dts); + activity=cat(1,activity,dactivity{NN,cue}(trls,:)); + end + + parfor t=1:tps + predictionO=NaN(dts*6,1); + predictionV=NaN(dts*6,1); + predictionS=NaN(dts*6,1); + for fold=1:folds + sel=trains(:,fold)==1; + aoi=activity(sel,t); + if sum(aoi)>0 %only run if there are spikes + mdl = fitcdiscr(aoi,odorLabels(sel)); + predictionO(sel==0)=predict(mdl,activity(sel==0,t)); + + mdl = fitcdiscr(aoi,valueLabels(sel)); + predictionV(sel==0)=predict(mdl,activity(sel==0,t)); + + mdl = fitcdiscr(aoi,setLabels(sel)); + predictionS(sel==0)=predict(mdl,activity(sel==0,t)); + end + end + difference=predictionO-odorLabels; + odorPerf(NN,t)=sum(difference==0)/length(difference); + + difference=predictionV-valueLabels; + valuePerf(NN,t)=sum(difference==0)/length(difference); + + difference=predictionS-setLabels; + setPerf(NN,t)=sum(difference==0)/length(difference); + end + if rem(NN,10)==0 fprintf('Neuron %d \n',NN); end +end +toc +save('dactivity.mat','dactivity','odorPerf','valuePerf','setPerf','-v7.3'); +end +%% plot single neuron decoding of cue identity +cats={1,2,3}; +catcolors={[0 0.7 1],[0.7 0 1],[0.6 0.6 0.6],[0.05 0.25 0.45]}; +vdata=[]; +vcat=[]; +vsess=[]; +vtp=[]; +odata=[]; +for c=1:length(cats) + sel=ismember(category,cats{c}); + + subplot(1,3,1) + hold on + up=nanmean(valuePerf(sel,:),1)+nanste(valuePerf(sel,:),1); + down=nanmean(valuePerf(sel,:),1)-nanste(valuePerf(sel,:),1); + plot(dbinTimes,nanmean(valuePerf(sel,:),1),'color',catcolors{c},'linewidth',1); + patch([dbinTimes,dbinTimes(end:-1:1)],[up,down(end:-1:1)],catcolors{c},'EdgeColor','none');alpha(0.5); + plot(dwindow,[1/3 1/3],':','color','k','linewidth',0.75) + + title('Single unit decoding of value'); + xlim(dwindow); + ylim([0.3 0.5]); + + ylabel('Decoding accuracy'); + xlabel('Time from odor onset (s)'); + + vdata=cat(1,vdata,valuePerf(sel,:)); + vcat=cat(1,vcat,c*ones(sum(sel),1)); + vsess=cat(1,vsess,neuronSession(sel)); + vtp=cat(1,vtp,repmat(1:size(valuePerf,2),sum(sel),1)); + + subplot(1,3,2) + hold on + up=nanmean(odorPerf(sel,:),1)+nanste(odorPerf(sel,:),1); + down=nanmean(odorPerf(sel,:),1)-nanste(odorPerf(sel,:),1); + plot(dbinTimes,nanmean(odorPerf(sel,:),1),'color',catcolors{c},'linewidth',1); + patch([dbinTimes,dbinTimes(end:-1:1)],[up,down(end:-1:1)],catcolors{c},'EdgeColor','none');alpha(0.5); + plot(dwindow,[1/6 1/6],':','color','k','linewidth',0.75) + + title('Single unit decoding of odor'); + xlim(dwindow); + ylim([0.15 0.27]); + + odata=cat(1,odata,odorPerf(sel,:)); + +end + +%value anova +vcat=repmat(vcat,1,size(valuePerf,2)); +vsess=repmat(vsess,1,size(valuePerf,2)); + +%odor anova +[~,tbl,stats]=anovan(odata(:),{categorical(vcat(:)),categorical(vtp(:)),categorical(vsess(:))},'varnames',{'cat','time','sess'},'random',[3],'display','off',... + 'model',[1 0 0;0 1 0;0 0 1;1 1 0]); +c=multcompare(stats,'dimension',[1 2],'ctype','bonferroni'); +%% population decoding of cue identity using different cue types +numN=[1 5 10 25 50 75 100 200]; +pwindow=[1 2.5]; +tpsel=dbinTimes>pwindow(1) & dbinTimes0 + nactivity=nactivity/max(nactivity,[],'all'); + end + activity(:,n)=nactivity; + end + + predictionO=NaN(dts*6,1); + predictionV=NaN(dts*6,1); + for fold=1:folds + sel=trains(:,fold)==1; + aoi=activity(sel,:); + aot=activity(sel==0,:); + if any(sum(aoi,1)==0) %only include neuron if there are spikes + aot=activity(sel==0,sum(aoi,1)>0); + aoi=activity(sel,sum(aoi,1)>0); + end + + mdl = fitcdiscr(aoi,odorLabels(sel),'gamma',gamma); + predictionO(sel==0)=predict(mdl,aot); + + mdl = fitcdiscr(aoi,valueLabels(sel),'gamma',gamma); + predictionV(sel==0)=predict(mdl,aot); + + end + difference=predictionO-odorLabels; + odorPerfP(r,nn,c)=sum(difference==0)/length(difference); + + difference=predictionV-valueLabels; + valuePerfP(r,nn,c)=sum(difference==0)/length(difference); + + + end + fprintf('Category %d, %d neurons \n',c,n2u); + end +end + +save('dactivity.mat','dactivity','odorPerf','valuePerf','setPerf','odorPerfP','valuePerfP','-v7.3'); + +end + +%% plot population decoding of cue identity (by neuron category) +figure; +cats={1,2,3}; +catcolors={[0 0.7 1],[0.7 0 1],[0.6 0.6 0.6],[0.05 0.25 0.45]}; + +subplot(1,3,1) +hold on +for c=1:length(cats) + errorbar([1:length(numN)]-0.3+0.15*c,nanmean(squeeze(valuePerfP(:,:,c)),1),nanstd(squeeze(valuePerfP(:,:,c)),1),'o','color',catcolors{c},'linewidth',1); +end +plot([0.5 length(numN)+0.5],[1/3 1/3],':','color','k','linewidth',0.75) +ylabel('Decoding accuracy'); +xlabel('Neurons used'); +ylim([0 1.1]); +xlim([0.5 length(numN)+0.5]); +yticks(0:0.25:1); +xticks(1:length(numN)); +xticklabels(numN); +title('Value'); + +%pvals +pvals={}; +pvalCh=NaN(length(cats)*2,length(numN)); +for n=1:length(numN) + for c1=1:length(cats) + for c2=1:length(cats) + pvals{1,n}(c1,c2)=(sum(squeeze(valuePerfP(:,n,c1))<=squeeze(valuePerfP(:,n,c2))','all')+1)/(size(valuePerfP,1)^2+1); + end + pvalCh(c1,n)=(sum(squeeze(valuePerfP(:,n,c1))<=(1/3))+1)/(size(valuePerfP,1)+1); + + end +end + +subplot(1,3,2) +hold on +for c=1:length(cats) + errorbar([1:length(numN)]-0.25+0.1*c,nanmean(squeeze(odorPerfP(:,:,c)),1),nanstd(squeeze(valuePerfP(:,:,c)),1),'o','color',catcolors{c},'linewidth',1); +end +plot([0.5 length(numN)+0.5],[1/6 1/6],':','color','k','linewidth',0.75) +ylim([0 1.1]); +xlim([0.5 length(numN)+0.5]); +yticks(0:0.25:1); +xticks(1:length(numN)); +xticklabels(numN); +title('Odor'); + +%pvals +for n=1:length(numN) + for c1=1:length(cats) + for c2=1:length(cats) + pvals{2,n}(c1,c2)=(sum(squeeze(odorPerfP(:,n,c1))<=squeeze(odorPerfP(:,n,c2))','all')+1)/(size(valuePerfP,1)^2+1); + end + pvalCh(c1+3,n)=(sum(squeeze(odorPerfP(:,n,c1))<=(1/6))+1)/(size(valuePerfP,1)+1); + + end +end + + +%% population decoding of value using different cue types +numN=[1 5 10 25 50 75 100 200]; +pwindow=[1 2.5]; +tpsel=dbinTimes>pwindow(1) & dbinTimes0 + end + activity(:,n)=nactivity; + end + + predictionV=NaN(dts*6,1); + for fold=1:folds + trainsel=trains(:,fold)==1; + testsel=tests(:,fold)==1; + aoi=activity(trainsel,:); + aot=activity(testsel,:); + if any(sum(aoi,1)==0) %only include neuron if there are spikes + aot=activity(testsel,sum(aoi,1)>0); + aoi=activity(trainsel,sum(aoi,1)>0); + end + + fit = glmnet(aoi,valueLabels(trainsel),'gaussian',glmopts); + prediction = glmnetPredict(fit,aot,[],'response'); + predictionV(testsel)=prediction(:,end); + + + end + + catPrediction=NaN(length(predictionV),1); + catPrediction(predictionV<0.25 & predictionV>=-0.25)=0; + catPrediction(predictionV<0.75 & predictionV>=0.25)=0.5; + catPrediction(predictionV<1.25 & predictionV>=0.75)=1; + + difference=catPrediction-valueLabels; + valuePerfP(r,nn,c)=sum(difference==0)/length(difference); + + for cue=1:3 + valueCueP(r,nn,c,cue)=mean(predictionV((cue-1)*100+1:cue*100)); + end + + end + fprintf('Category %d, %d neurons \n',c,n2u); + end +end + +save('dactivityV.mat','dactivity','valuePerf','valuePerfP','-v7.3'); + +end + +%% plot population decoding of value (by neuron category) +figure; +cats={1,2,3}; +catcolors={[0 0.7 1],[0.7 0 1],[0.6 0.6 0.6],[0.05 0.25 0.45]}; + +subplot(1,3,1) +hold on +for c=1:length(cats) + errorbar([1:length(numN)]-0.3+0.15*c,nanmean(squeeze(valuePerfP(:,:,c)),1),nanstd(squeeze(valuePerfP(:,:,c)),1),'o','color',catcolors{c},'linewidth',1); +end +plot([0.5 length(numN)+0.5],[1/3 1/3],':','color','k','linewidth',0.75) +ylabel('Decoding accuracy'); +xlabel('Neurons used'); +ylim([0 1.1]); +xlim([0.5 length(numN)+0.5]); +yticks(-1:0.5:1); +xticks(1:length(numN)); +xticklabels(numN); +title('Value'); + + subplot(2,3,2) + hold on +for c=1:length(cats) + errorbar([1:length(numN)]-0.3+0.15*c,nanmean(squeeze(valueCueP(:,:,c,1)),1),nanstd(squeeze(valueCueP(:,:,c,1)),1),'o','color',catcolors{c},'linewidth',1); +end +plot([0.5 length(numN)+0.5],[1 1],':','color',colors{1},'linewidth',0.75) +plot([0.5 length(numN)+0.5],[0.5 0.5],':','color',colors{2},'linewidth',0.75) + +ylabel('Predicted cue value'); +ylim([0.4 1.1]); +xlim([0.5 length(numN)+0.5]); +yticks(-1:0.5:1); +xticks([]); +xticklabels(numN); +title('CS+'); + + subplot(2,3,5) + hold on +for c=1:length(cats) + errorbar([1:length(numN)]-0.3+0.15*c,nanmean(squeeze(valueCueP(:,:,c,3)),1),nanstd(squeeze(valueCueP(:,:,c,3)),1),'o','color',catcolors{c},'linewidth',1); +end +plot([0.5 length(numN)+0.5],[0 0],':','color',colors{3},'linewidth',0.75) +plot([0.5 length(numN)+0.5],[0.5 0.5],':','color',colors{2},'linewidth',0.75) + +xlabel('Neurons used'); +ylim([-0.1 0.6]); +xlim([0.5 length(numN)+0.5]); +yticks(-1:0.5:1); +xticks(1:length(numN)); +xticklabels(numN); +title('CS-'); + +%pvals +pvals={}; +pvalCh=NaN(length(cats)*2,length(numN)); +for n=1:length(numN) + for c1=1:length(cats) + for c2=1:length(cats) + pvals{1,n}(c1,c2)=(sum(squeeze(valuePerfP(:,n,c1))<=squeeze(valuePerfP(:,n,c2))','all')+1)/(size(valuePerfP,1)^2+1); + end + pvalCh(c1,n)=(sum(squeeze(valuePerfP(:,n,c1))<=(1/3))+1)/(size(valuePerfP,1)+1); + + end +end + +%% population decoding of value using different regions +n2u=5; +pwindows={[-0.5 0],[1 2.5]}; +reps=1000; +preRun=false; + +regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; +neuronRegionOlf=neuronRegionAdj; +neuronRegionOlf(ismember(neuronRegionAdj,{'EPd','PIR'}))={'OLF'}; + +if preRun==false + + +folds=6; +tests=zeros(dts*6,folds); +for f=1:folds +tests((f-1)*dts+1:f*dts,f)=1; +end + +trains1=zeros(dts*6,1);trains1(51:100)=1;trains1(126:175)=1;trains1(226:275)=1; +trains2=zeros(dts*6,1);trains2(1:50)=1;trains2(126:175)=1;trains2(226:275)=1; +trains3=zeros(dts*6,1);trains3(26:75)=1;trains3(151:200)=1;trains3(226:275)=1; +trains4=zeros(dts*6,1);trains4(26:75)=1;trains4(101:150)=1;trains4(226:275)=1; +trains5=zeros(dts*6,1);trains5(26:75)=1;trains5(126:175)=1;trains5(251:300)=1; +trains6=zeros(dts*6,1);trains6(26:75)=1;trains6(126:175)=1;trains6(201:250)=1; +trains=[trains1 trains2 trains3 trains4 trains5 trains6]; + +valueLabels=[ones(dts,1);ones(dts,1);0.5*ones(dts,1);0.5*ones(dts,1);zeros(dts,1);zeros(dts,1)]; + +tps=size(dactivity{1,1},2); + +valuePerfR=NaN(reps,length(regions),length(pwindows)); + +for tp=1:length(pwindows) + tpsel=dbinTimes>pwindows{tp}(1) & dbinTimes0 + nactivity=nactivity/max(nactivity,[],'all'); + end + activity(:,n)=nactivity; + end + + predictionV=NaN(dts*6,1); + for fold=1:folds + trainsel=trains(:,fold)==1; + testsel=tests(:,fold)==1; + aoi=activity(trainsel,:); + aot=activity(testsel,:); + if any(sum(aoi,1)==0) %only include neuron if there are spikes + aot=activity(testsel,sum(aoi,1)>0); + aoi=activity(trainsel,sum(aoi,1)>0); + end + + fit = glmnet(aoi,valueLabels(trainsel),'gaussian',glmopts); + prediction = glmnetPredict(fit,aot,[],'response'); + predictionV(testsel)=prediction(:,end); + + end + + catPrediction=NaN(length(predictionV),1); + catPrediction(predictionV<0.25 & predictionV>=-0.25)=0; + catPrediction(predictionV<0.75 & predictionV>=0.25)=0.5; + catPrediction(predictionV<1.25 & predictionV>=0.75)=1; + + difference=catPrediction-valueLabels; + valuePerfR(r,c,tp)=sum(difference==0)/length(difference); + + end + fprintf('TP %d, region %d \n',tp,c); + end +end +save('dactivityV.mat','dactivity','valuePerf','valuePerfP','valuePerfR','-v7.3'); + +end +%% plot population decoding of value (by region) +figure; +catcolors={[0 0 0],[0.35 0.35 1]}; +catcolors={[0 0 0],[0 0.7 1]}; + +subplot(1,3,1) +hold on +for tp=1:length(pwindows) +errorbar(1:length(regions),nanmean(squeeze(valuePerfR(:,:,tp)),1),nanstd(squeeze(valuePerfR(:,:,tp)),1),'o','color',catcolors{tp},'linewidth',1); +end +plot([0.5 length(regions)+0.5],[1/3 1/3],':','color','k','linewidth',0.75) +ylabel('Decoding accuracy'); +xlabel('Regions'); +ylim([0 1.1]); +xlim([0.5 length(regions)+0.5]); +yticks(0:0.25:1); +xticks(1:length(regions)); +xticklabels(regions); +xtickangle(45); +title('Value'); + +%pvals +pvalsV=[]; +pvalB=NaN(2,length(regions)); +for c1=1:length(regions) + for c2=1:length(regions) + pvalsV(c1,c2)=(sum(squeeze(valuePerfR(:,c1,2))<=squeeze(valuePerfR(:,c2,2))','all')+1)/(size(valuePerfR,1)^2+1); + end + pvalsB(1,c1)=(sum(squeeze(valuePerfR(:,c1,2))<=squeeze(valuePerfR(:,c1,1))','all')+1)/(size(valuePerfR,1)^2+1); + + +end + + +%% population decoding of value using different region groups +n2u=25; +pwindows={[-0.5 0],[1 2.5]}; +reps=1000; +preRun=false; + +reggrps=[1:3]; +if preRun==false + +dts=50; +gamma=1; + +folds=6; +tests=zeros(dts*6,folds); +for f=1:folds +tests((f-1)*dts+1:f*dts,f)=1; +end + +trains1=zeros(dts*6,1);trains1(51:100)=1;trains1(126:175)=1;trains1(226:275)=1; +trains2=zeros(dts*6,1);trains2(1:50)=1;trains2(126:175)=1;trains2(226:275)=1; +trains3=zeros(dts*6,1);trains3(26:75)=1;trains3(151:200)=1;trains3(226:275)=1; +trains4=zeros(dts*6,1);trains4(26:75)=1;trains4(101:150)=1;trains4(226:275)=1; +trains5=zeros(dts*6,1);trains5(26:75)=1;trains5(126:175)=1;trains5(251:300)=1; +trains6=zeros(dts*6,1);trains6(26:75)=1;trains6(126:175)=1;trains6(201:250)=1; +trains=[trains1 trains2 trains3 trains4 trains5 trains6]; + +valueLabels=[ones(dts,1);ones(dts,1);0.5*ones(dts,1);0.5*ones(dts,1);zeros(dts,1);zeros(dts,1)]; + +valuePerfG=NaN(reps,length(reggrps),length(pwindows)); + +for tp=1:length(pwindows) + tpsel=dbinTimes>pwindows{tp}(1) & dbinTimes0 + nactivity=nactivity/max(nactivity,[],'all'); + end + activity(:,n)=nactivity; + end + + predictionV=NaN(dts*6,1); + for fold=1:folds + trainsel=trains(:,fold)==1; + testsel=tests(:,fold)==1; + aoi=activity(trainsel,:); + aot=activity(testsel,:); + if any(sum(aoi,1)==0) %only include neuron if there are spikes + aot=activity(testsel,sum(aoi,1)>0); + aoi=activity(trainsel,sum(aoi,1)>0); + end + + fit = glmnet(aoi,valueLabels(trainsel),'gaussian',glmopts); + prediction = glmnetPredict(fit,aot,[],'response'); + predictionV(testsel)=prediction(:,end); + + end + + catPrediction=NaN(length(predictionV),1); + catPrediction(predictionV<0.25 & predictionV>=-0.25)=0; + catPrediction(predictionV<0.75 & predictionV>=0.25)=0.5; + catPrediction(predictionV<1.25 & predictionV>=0.75)=1; + + difference=catPrediction-valueLabels; + valuePerfG(r,c,tp)=sum(difference==0)/length(difference); + + end + fprintf('TP %d, region %d \n',tp,c); + end +end +save('dactivityV.mat','dactivity','valuePerf','valuePerfP','valuePerfR','valuePerfG','-v7.3'); + +end +%% plot population decoding of value (by region group) +figure; +catcolors={[0 0 0],[0.35 0.35 1]}; +catcolors={[0 0 0],[0 0.7 1]}; + +subplot(1,6,1) +hold on +for tp=1:length(pwindows) +errorbar(1:length(reggrps),nanmean(squeeze(valuePerfG(:,:,tp)),1),nanstd(squeeze(valuePerfG(:,:,tp)),1),'o','color',catcolors{tp},'linewidth',1); +end +plot([0.5 length(reggrps)+0.5],[1/3 1/3],':','color','k','linewidth',0.75) +ylabel('Decoding accuracy'); +xlabel('Region groups'); +ylim([0 1.1]); +xlim([0.5 length(reggrps)+0.5]); +yticks(0:0.25:1); +xticks(1:length(reggrps)); +xticklabels(regionGroupNames(reggrps)); +xtickangle(45); +title('Value'); + +%pvals +pvalsV=[]; +for c1=1:length(reggrps) + for c2=1:length(reggrps) + pvalsV(c1,c2)=(sum(squeeze(valuePerfG(:,c1,2))<=squeeze(valuePerfG(:,c2,2))','all')+1)/(size(valuePerfG,1)^2+1); + end +end + diff --git a/Code/imagingAcquisition.m b/Code/imagingAcquisition.m index 3c88246..4251802 100644 --- a/Code/imagingAcquisition.m +++ b/Code/imagingAcquisition.m @@ -1148,6 +1148,122 @@ yticks([]); end xlabel('seconds from odor onset'); +%% correlating activity from days 1 to 2 to 3 + +plotBins=binTimes>=analysisWindow(1) & binTimes<=analysisWindow(2); + +map=redblue(256); +figure; +colormap(map); + +activity1 = cat(2,dffPSTH3.o1d1{:,3}); +activity1 = activity1(trackedROIids{1}>0,plotBins); +[~,idOrder]=sort(trackedROIids{1}(trackedROIids{1}>0)); +activity1=activity1(idOrder,:); + +activity2 = cat(2,dffPSTH3.o1d2{:,3}); +activity2 = activity2(trackedROIids{2}>0,plotBins); +[~,idOrder]=sort(trackedROIids{2}(trackedROIids{2}>0)); +activity2=activity2(idOrder,:); + +activity3 = cat(2,dffPSTH3.o1d3{:,3}); +activity3 = activity3(trackedROIids{3}>0,plotBins); +[~,idOrder]=sort(trackedROIids{3}(trackedROIids{3}>0)); +activity3=activity3(idOrder,:); + +numneur=size(activity1,1); +[~,randomize]=sort(activity2(:,2)); + +corr12=corrcoef([activity1' activity2']); +corr12=corr12(numneur+1:end,1:numneur); +corr12sh=corrcoef([activity1' activity2(randomize,:)']); +corr12sh=corr12sh(numneur+1:end,1:numneur); +percentile=NaN(numneur,1); +percentileSh=NaN(numneur,1); +for n=1:numneur + self=corr12(n,n); + other=corr12(n,:);other(n)=[]; + percentile(n,1)=sum(other<=self)/numneur; + + self=corr12sh(n,n); + other=corr12sh(n,:);other(n)=[]; + percentileSh(n,1)=sum(other<=self)/numneur; +end + +subplot(1,3,1) +hold on +[y,x]=ecdf(percentileSh*100); +plot(x,y,'color',[0.6 0.6 0.6],'linewidth',1); +[y,x]=ecdf(percentile*100); +plot(x,y,'color','k','linewidth',1); +xlabel('Percentile'); +ylabel('Cumulative fraction'); +title('Correlation rank'); +legend('Shuffle','True'); + + +subplot(1,3,3) +hold on +mouseRank=NaN(length(mice),2); +for mouse=1:length(mice) + mouseSel=strcmp(trackedMouse,mice(mouse)); + mouseRank(mouse,1)=median(percentile(mouseSel)); + mouseRank(mouse,2)=median(percentileSh(mouseSel)); +end + +scatter(rand(length(mice),1)/2+0.5,mouseRank(:,1),24,'k','filled'); +scatter(rand(length(mice),1)/2+0.5,mouseRank(:,2),24,[0.6 0.6 0.6],'filled'); +errorbar(1.25,mean(mouseRank(:,1)),nanste(mouseRank(:,1),1),'o','color','k','linewidth',1); +errorbar(1.25,mean(mouseRank(:,2)),nanste(mouseRank(:,2),1),'o','color',[0.6 0.6 0.6],'linewidth',1); + +corr23=corrcoef([activity2' activity3']); +corr23=corr23(numneur+1:end,1:numneur); +corr23sh=corrcoef([activity2(randomize,:)' activity3']); +corr23sh=corr23sh(numneur+1:end,1:numneur); + +percentile=NaN(numneur,1); +for n=1:numneur + self=corr23(n,n); + other=corr23(n,:);other(n)=[]; + percentile(n,1)=sum(other<=self)/numneur; + + self=corr23sh(n,n); + other=corr23sh(n,:);other(n)=[]; + percentileSh(n,1)=sum(other<=self)/numneur; +end + +subplot(1,3,2) +hold on +[y,x]=ecdf(percentileSh*100); +plot(x,y,'color',[0.6 0.6 0.6],'linewidth',1); +[y,x]=ecdf(percentile*100); +plot(x,y,'color','k','linewidth',1); +xlabel('Percentile'); +ylabel('Cumulative fraction'); +title('Correlation rank'); +legend('shuffled','real'); + +subplot(1,3,3) +hold on +mouseRank=NaN(length(mice),2); +for mouse=1:length(mice) + mouseSel=strcmp(trackedMouse,mice(mouse)); + mouseRank(mouse,1)=median(percentile(mouseSel)); + mouseRank(mouse,2)=median(percentileSh(mouseSel)); +end + +scatter(rand(length(mice),1)/2+2.5,mouseRank(:,1),24,'k','filled'); +scatter(rand(length(mice),1)/2+2.5,mouseRank(:,2),24,[0.6 0.6 0.6],'filled'); +errorbar(3.25,mean(mouseRank(:,1)),nanste(mouseRank(:,1),1),'o','color','k','linewidth',1); +errorbar(3.25,mean(mouseRank(:,2)),nanste(mouseRank(:,2),1),'o','color',[0.6 0.6 0.6],'linewidth',1); +axis([0 4 0 1]); +xticks([1 3]); +yticks([0 0.5 1]); +xticklabels({'day 1->2','day 2->3'}); +ylabel('Correlation rank'); + + + %% glm performance across tracked cells %starter model @@ -1318,6 +1434,200 @@ end +%% plotting activity of CS+ preferring neurons day 1 to 2 to 3 + +%sort based on category on day 3 +catInc=category{3}(trackedROIids{3}>0); +[~,idOrder]=sort(trackedROIids{3}(trackedROIids{3}>0)); +catTrack=catInc(idOrder); + +plotWindow=[-0.5 2.5]; +plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); +xvals=binTimes(plotBins); + +%also select peak activity for CS+ +cueWindow = [0 2.5]; +cueBins=binTimes>=cueWindow(1) & binTimes<=cueWindow(2); +cspActivity=mean(dffPSTH3.o1d3{1,3}((trackedROIids{3}>0),cueBins),2); +cueRespP=abs(mean(dffPSTH3.o1d3{1,3}(:,cueBins),2)); +cueRespF=abs(mean(dffPSTH3.o1d3{2,3}(:,cueBins),2)); +cueRespM=abs(mean(dffPSTH3.o1d3{3,3}(:,cueBins),2)); +CSPpref = cueRespP > cueRespF & cueRespP > cueRespM; +cueInc = CSPpref(trackedROIids{3}>0); +cueTrack = cueInc(idOrder); +sel = catTrack==1 & cueTrack; + +sessionStage=3; +colors{1,1}=[0.1 0.6 0.2]; +colors{2,1}=[0.4 0.1 0.4]; +colors{3,1}=[0.3 0.3 0.3]; +directions=[1 -1]; +specs={'-','--'}; +diractivity = dffPSTH6.o1d3{(1-1)*2+2,1}(:,plotBins); +direction=sign(max(diractivity,[],2)-abs(min(diractivity,[],2))); +direction=direction(trackedROIids{3}>0); +direction=direction(idOrder); +for session=1:length(sessions) + sn=sessions{session}; + for d=1:2 + subplot(2,3,(d-1)*3+session); + + hold on; + for cue=1:3 + + if session==3 activity = dffPSTH6.(sn){(cue-1)*2+1}(trackedROIids{session}>0,:); + else activity = dffPSTH3.(sn){cue,sessionStage}(trackedROIids{session}>0,:); end + + [~,idOrder]=sort(trackedROIids{session}(trackedROIids{session}>0)); + activity=activity(idOrder,:); + + %get values + psth=nanmean(activity(sel&direction==directions(d),plotBins)); + sem=nanste(activity(sel&direction==directions(d),plotBins),1); %calculate standard error of the mean + up=psth+sem; + down=psth-sem; + + %plotting + + plot(binTimes(plotBins),psth,'Color',colors{cue,1},'linewidth',0.75); + patch([xvals,xvals(end:-1:1)],[up,down(end:-1:1)],colors{cue,1},'EdgeColor','none');alpha(0.2); + + + end + + + + if d==1 text(0.1,1.5,num2str(sum(sel&direction==directions(d)))); end + if d==2 text(0.1,-0.45,num2str(sum(sel&direction==directions(d)))); end + if d==2 + xlabel('seconds from odor onset'); + ylabel('z-score'); + else + xticks([]); + end + + %if reg>1 yticks([]); end + plot([0 0],[-1 5],'color','k','linewidth',0.75); + plot(plotWindow,[0 0],':','color','k','linewidth',0.5); + if d==1 axis([plotWindow -0.25 4]); end + if d==2 axis([plotWindow -0.8 0.4]); yticks([-0.5 0]); end + end +end + + +%% correlating activity from days 1 to 2 to 3 for CS+preferring cue neurons (run immediately after code above) + +plotBins=binTimes>=analysisWindow(1) & binTimes<=analysisWindow(2); + +map=redblue(256); +figure; +colormap(map); + +activity1 = cat(2,dffPSTH3.o1d1{:,3}); +activity1 = activity1(trackedROIids{1}>0,plotBins); +[~,idOrder]=sort(trackedROIids{1}(trackedROIids{1}>0)); +activity1=activity1(idOrder,:); + +activity2 = cat(2,dffPSTH3.o1d2{:,3}); +activity2 = activity2(trackedROIids{2}>0,plotBins); +[~,idOrder]=sort(trackedROIids{2}(trackedROIids{2}>0)); +activity2=activity2(idOrder,:); + +activity3 = cat(2,dffPSTH3.o1d3{:,3}); +activity3 = activity3(trackedROIids{3}>0,plotBins); +[~,idOrder]=sort(trackedROIids{3}(trackedROIids{3}>0)); +activity3=activity3(idOrder,:); + +numneur=size(activity1,1); +[~,randomize]=sort(activity2(:,2)); + +corr12=corrcoef([activity1' activity2']); +corr12=corr12(numneur+1:end,1:numneur); +corr12sh=corrcoef([activity1' activity2(randomize,:)']); +corr12sh=corr12sh(numneur+1:end,1:numneur); +percentile=NaN(numneur,1); +percentileSh=NaN(numneur,1); +for n=1:numneur + self=corr12(n,n); + other=corr12(n,:);other(n)=[]; + percentile(n,1)=sum(other<=self)/numneur; + + self=corr12sh(n,n); + other=corr12sh(n,:);other(n)=[]; + percentileSh(n,1)=sum(other<=self)/numneur; +end + +subplot(1,3,1) +hold on +[y,x]=ecdf(percentileSh(sel)*100); +plot(x,y,'color',[0.6 0.6 0.6],'linewidth',1); +[y,x]=ecdf(percentile(sel)*100); +plot(x,y,'color','k','linewidth',1); +xlabel('Percentile'); +ylabel('Cumulative fraction'); +title('Correlation rank'); +legend('Shuffle','True'); + +subplot(1,3,3) +hold on +mouseRank=NaN(length(mice),2); +for mouse=1:length(mice) + mouseSel=strcmp(trackedMouse,mice(mouse))&sel; + mouseRank(mouse,1)=median(percentile(mouseSel)); + mouseRank(mouse,2)=median(percentileSh(mouseSel)); +end + +scatter(rand(length(mice),1)/2+0.5,mouseRank(:,1),24,'k','filled'); +scatter(rand(length(mice),1)/2+0.5,mouseRank(:,2),24,[0.6 0.6 0.6],'filled'); +errorbar(1.25,nanmean(mouseRank(:,1)),nanste(mouseRank(:,1),1),'o','color','k','linewidth',1); +errorbar(1.25,nanmean(mouseRank(:,2)),nanste(mouseRank(:,2),1),'o','color',[0.6 0.6 0.6],'linewidth',1); + +corr23=corrcoef([activity2' activity3']); +corr23=corr23(numneur+1:end,1:numneur); +corr23sh=corrcoef([activity2(randomize,:)' activity3']); +corr23sh=corr23sh(numneur+1:end,1:numneur); + +percentile=NaN(numneur,1); +for n=1:numneur + self=corr23(n,n); + other=corr23(n,:);other(n)=[]; + percentile(n,1)=sum(other<=self)/numneur; + + self=corr23sh(n,n); + other=corr23sh(n,:);other(n)=[]; + percentileSh(n,1)=sum(other<=self)/numneur; +end + +subplot(1,3,2) +hold on +[y,x]=ecdf(percentileSh(sel)*100); +plot(x,y,'color',[0.6 0.6 0.6],'linewidth',1); +[y,x]=ecdf(percentile(sel)*100); +plot(x,y,'color','k','linewidth',1); +xlabel('Percentile'); +ylabel('Cumulative fraction'); +title('Correlation rank'); +legend('shuffled','real'); + +subplot(1,3,3) +hold on +mouseRank=NaN(length(mice),2); +for mouse=1:length(mice) + mouseSel=strcmp(trackedMouse,mice(mouse))&sel; + mouseRank(mouse,1)=median(percentile(mouseSel)); + mouseRank(mouse,2)=median(percentileSh(mouseSel)); +end + +scatter(rand(length(mice),1)/2+2.5,mouseRank(:,1),24,'k','filled'); +scatter(rand(length(mice),1)/2+2.5,mouseRank(:,2),24,[0.6 0.6 0.6],'filled'); +errorbar(3.25,nanmean(mouseRank(:,1)),nanste(mouseRank(:,1),1),'o','color','k','linewidth',1); +errorbar(3.25,nanmean(mouseRank(:,2)),nanste(mouseRank(:,2),1),'o','color',[0.6 0.6 0.6],'linewidth',1); +axis([0 4 0 1]); +xticks([1 3]); +yticks([0 0.5 1]); +xticklabels({'day 1->2','day 2->3'}); +ylabel('Correlation rank'); + %% GLM task variables in spatial location diff --git a/Code/imagingOdorSetsValue.m b/Code/imagingOdorSetsValue.m index e6a6277..34e204f 100644 --- a/Code/imagingOdorSetsValue.m +++ b/Code/imagingOdorSetsValue.m @@ -691,23 +691,6 @@ lambda=[0.001 0.005 0.01 0.02 0.05 0.1 0.2 0.5]; binsPerCue=diff(windows{1})/binSize; -%permute all possible combinations of odors for meaning shuffle -allodors=1:6; -perms6=NaN(90,6); -first2=nchoosek(allodors,2); -perm=0; -for f2=1:length(first2) - remaining=allodors(~ismember(allodors,first2(f2,:))); - next2=nchoosek(remaining,2); - for n2=1:length(next2) - last=remaining(~ismember(remaining,next2(n2,:))); - perm=perm+1; - perms6(perm,:)=[first2(f2,:) next2(n2,:) last]; - end -end -perms6all=perms6; -varExpValueSh=NaN(totalNeurons,4+length(perms6all)); - NN=0; for mouse=1:length(mice) session1=mouse; @@ -727,43 +710,7 @@ As{sub} = At; rAs{sub} = At * bR(:,1:components); end - - %get original cue order for shuffle - firstBins=[]; - cueOrder=[1 3 5 2 4 6]; - for cue=1:6 - firstBins{cue}=find(A(:,1+(cue-1)*binsPerCue)); - firstBins{cue}=[firstBins{cue} ones(length(firstBins{cue}),1)*cueOrder(cue)]; - end - alltrials=cat(1,firstBins{:}); - [sortedtrials,ordered]=sort(alltrials(:,1)); - originalOrder=alltrials(ordered,2); - combinedValues=[trialValues{mouse,1};trialValues{mouse,2}]; - - %make new shuffled predictor matrices - A3Sh={}; - for perm=1:length(perms6all)+1 - if perm>size(perms6all,1) - newValue=ones(length(alltrials),1); - else - newOrder=NaN(length(alltrials),1); - for cue=1:6 - newOrder(originalOrder==cue)=perms6all(perm,cue); - end - - newValue=NaN(length(alltrials),1); - for cue=1:6 - originalValues=combinedValues(originalOrder==cue); - newValue(newOrder==cue)=mean(originalValues(randsample(sum(originalOrder==cue),sum(newOrder==cue),'true'))); - end - end - - A3Sh{perm}=XallC3{mouse}; - for bin=1:binsPerCue - A3Sh{perm}(sortedtrials(:,1)-1+bin,bin)=newValue; - end - end - + for roi=1:size(Yall{session1},2) NN=NN+1; y1=Yall{session1}(:,roi); @@ -775,9 +722,7 @@ for fold=1:folds train=trains{fold}; test=train==0; - % fit = glmnet(rA(train,:), y(train), 'gaussian', opts); - % this_a = glmnetCoef(fit, lambda); - % fitK = this_a(2:end); + fitK=lassoglm(rA(train,:),y(train),'normal','alpha',0.5,'lambda',lambda); kernel=bR(:,1:components)*fitK; predLam(test,:)=A(test,:)*kernel; @@ -812,54 +757,12 @@ end predF{NN,1+sub}=pred; varExp(NN,1+sub) = 1- var(y-pred)/var(y); - -% %variance over time -% for tb=1:binsPerTrial -% tbsel=tb:binsPerTrial:length(y); -% varExpT.(sn){NN,1+sub}(tb)=1- var(y(tbsel)-pred(tbsel))/var(y(tbsel)); -% end - end - - - - %models with alternate cue - %value instead of cue identitiy - A3 = XallC3{mouse}; - rA3 = A3 * bR3(:,1:valueComponents); - pred=NaN(size(y,1),1); - for fold=1:folds - train=trains{fold}; - test=train==0; - fitK=lassoglm(rA3(train,:),y(train),'normal','alpha',0.5,'lambda',thisLambda); - kernel=bR3(:,1:valueComponents)*fitK; - pred(test)=A3(test,:)*kernel; + end - varExp(NN,3+sub) = 1- var(y-pred)/var(y); - if varExp(NN,1)-varExp(NN,2)>0.02 - - for perm=1:size(perms6all,1)+1 - - rA3 = A3Sh{perm} * bR3(:,1:valueComponents); - trains = getfolds(rA3,folds,binsPerTrial); - - % if ismember(perm,topModels) - % display(perms6all(perm,:));display([combinedValues(1:10) newValue(1:10)]); - % end - - - pred=NaN(size(y,1),1); - for fold=1:folds - train=trains{fold}; - test=train==0; - fitK=lassoglm(rA3(train,:),y(train),'normal','alpha',0.5,'lambda',thisLambda); - kernel=bR3(:,1:valueComponents)*fitK; - pred(test)=A3Sh{perm}(test,:)*kernel; - end - varExpValueSh(NN,3+perm) = 1- var(y-pred)/var(y); - end + - end + end @@ -868,14 +771,78 @@ end -varExpValueSh(:,1)=varExp(:,1); -varExpValueSh(:,2)=varExp(:,2); -varExpValueSh(:,3)=varExp(:,5); +%% perform value analysis with a ton of shuffles +runglm=true; +%generate all value shuffles +valueShuffles = []; +valueVersion = []; +valueLevels = [1 1 0.5 0.5 0 0;1 1 1 1 0 0;1 1 1 0 0 0;1 1 0 0 0 0;1 0 0 0 0 0;1 1 1 1 1 0;1 1 1 1 1 1]; +for vl=1:size(valueLevels,1) +allCombos = perms(valueLevels(vl,:)); +allCombos = unique(allCombos,'rows'); +valueShuffles = cat(1,valueShuffles,allCombos); +valueVersion = cat(1,valueVersion,vl*ones(size(allCombos,1),1)); +end +valueShuffles=fliplr(valueShuffles); +folds=4; +NNstart=0; +varExpValueNew=NaN(totalNeurons,length(valueShuffles)); +%kernels=NaN(totalNeurons,31); +binsPerCue=diff(windows{1})/binSize; +opts.alpha=0.5; +glmopts=glmnetSet(opts); +if runglm +cueOrdering=[1 3 5 2 4 6]; +tic +for NS=1:length(mice) + + cueStarts={}; + for cue=1:6 + cueStarts{cue,1}(:,1)=find(XallC{NS}(:,(cue-1)*binsPerCue+1)); + cueStarts{cue,1}(:,2)=cueOrdering(cue)*ones(length(cueStarts{cue,1}),1); + end + + sortedStarts=sortrows(cat(1,cueStarts{:})); + cueOrder=sortedStarts(:,2); + + A = XallC{NS}(:,[1:binsPerCue size(XallC{NS},2)-1:size(XallC{NS},2)]); + trains = getfolds(A,folds,binsPerTrial); + for s=1:length(valueShuffles) + for c=1:length(cueOrder) + cueValue=valueShuffles(s,cueOrder(c)); + dm=cueValue*eye(binsPerCue); + A(c*binsPerTrial-binsPerCue+1:c*binsPerTrial,1:binsPerCue)=dm; + end + NN=NNstart; + for neuron=1:size(Yall{NS},2) + NN=NN+1; + y=YallC{NS}(:,neuron); + %cross-validated variance explained + pred=NaN(size(y,1),1); + for fold=1:folds + train=trains{fold}; + test=train==0; + %kernel=lassoglm(A(train,:),y(train),'normal','alpha',0.5,'lambda',thisLambda); + fit = glmnet(A(train,:),y(train),[],glmopts); + prediction = glmnetPredict(fit,A(test,:)); + pred(test)=prediction(:,end); + end + varExpValueNew(NN,s) = 1- var(y-pred)/var(y); +% kernel = fit.beta(:,end); +% if s==1 kernels(NN,:)=kernel; end + end + end + NNstart=NN; + fprintf('Mouse #%d \n',NS); +end +toc + +end -%% value coding analysis 90 with untuned +%% value coding analysis improvCutoff=0.02; overallCutoff=0.02; improvement=varExp(:,1)-varExp; @@ -886,123 +853,131 @@ cueK = improvement(:,2)>improvCutoff; %if best model is at least cutoff better than no cue kernels totalCueVariance=varExp(:,1)-varExp(:,2); %best model compared to no cue model -valueImp=varExp(:,1)-varExpValueSh; %how much better the best model is than value, shuffles, and one kernel +valueImp=varExp(:,1)-varExpValueNew; %how much better the best model is than value, shuffles, and one kernel %correct it so you can't be worse than no cues at all (should I do this?) for n=1:length(valueImp) valueImp(n,valueImp(n,:)>totalCueVariance(n))=totalCueVariance(n); end -%value specific variance -%whatever is accounted for by value kernel -cueVariance=totalCueVariance - valueImp(:,[4:end]); -[~,bestModel]=max(cueVariance,[],2); -bestModel(isnan(cueVariance(:,1)))=NaN; +includedModels = 1:length(valueShuffles); +%includedModels = [1:90 181:length(valueShuffles)]; +%includedModels = [1:length(valueShuffles)]; +shuffleLabels = cell(length(includedModels),1); +for sl=1:length(includedModels) + m=includedModels(sl); + shuffleLabels{sl} = sprintf('%g,%g,%g,%g,%g,%g',valueShuffles(m,1),valueShuffles(m,2),valueShuffles(m,3),valueShuffles(m,4),valueShuffles(m,5),valueShuffles(m,6)); +end +[~,bestModel]=max(varExpValueNew(:,includedModels),[],2); +% [~,bestModel]=max(varExpValueNew,[],2); +% bestModel(bestModel>90)=bestModel(bestModel>90)-90; +totalModels = length(includedModels); figure; -cueNames={'+','+','50','50','-','-'}; -odorValues=[0.5 0.5 0.37 0.37 0.05 0.05]; -vnames={}; -shuffVals=NaN(90,6); -for n=1:length(perms6all) - vnames{n,1}=append(cueNames{perms6all(n,:)==1},'/',cueNames{perms6all(n,:)==2},', ',cueNames{perms6all(n,:)==3},'/',cueNames{perms6all(n,:)==4},', ',cueNames{perms6all(n,:)==5},'/',cueNames{perms6all(n,:)==6}); - for odor=1:6 - shuffVals(n,odor)=odorValues(perms6all(n,:)==odor); - end -end +corrWVal=corrcoef(valueShuffles(includedModels,:)'); +corrWVal(isnan(corrWVal))=0; +[valCorr,valRank] = sort(corrWVal(:,1),'descend'); +[~,mdlPosition] = sort(valRank); -subplot(3,1,1); +subplot(5,1,1); hold on -frequency=histcounts(bestModel(cueK&~behavior&~isnan(bestModel)),0.5:1:91.5,'normalization','probability'); -outliers=frequency>0.05; -topModels=find(outliers); -b=bar([1:max(bestModel)-1 max(bestModel)+1],frequency,'facecolor',[0 0.2 0.4],'linewidth',0.01); +frequency=histcounts(bestModel(cueK&~behavior&~isnan(bestModel)),0.5:1:totalModels+0.5,'normalization','probability'); +%b=bar(1:max(bestModel),frequency,'facecolor',[0 0.2 0.4],'linewidth',0.01); +b=bar(mdlPosition,frequency,'facecolor',[0 0.2 0.4],'linewidth',0.01); + b.FaceColor='flat'; -for e=1:length(topModels) -b.CData(topModels(e),:)=[0.7 0 1]; %trial type +for e=2:17 +b.CData(e,:)=[0.7 0 1]; %trial type end b.CData(1,:)=[0 0.7 1]; %value -b.CData(end,:)=[0.6 0.6 0.6]; %non-specific - +b.CData(mdlPosition(end),:)=[0.6 0.6 0.6]; %non-specific +ylabel('fraction of cue cells'); +chance=1/totalModels; +plot([0 totalModels+1],[chance chance],':','color','k','linewidth',0.5); +xlim([0 totalModels+1]); +ylim([0 0.2]); +%xticks(mdlPosition(topModels(1:end-1))); +xticks([]); +yticks(0:0.1:0.3); +xlabel('cue coding models, sorted by similarity to value model'); +subplot(4,5,6); +hold on +scatter(valCorr(2:17,1),frequency(valRank(2:17)),24,[0.7 0 1],'filled'); +scatter(valCorr(1,1),frequency(valRank(1)),24,[0 0.7 1],'filled'); +scatter(valCorr(18:end,1),frequency(valRank(18:end)),24,[0 0.2 0.4],'filled'); +scatter(valCorr(mdlPosition(end),1),frequency(end),24,[0.6 0.6 0.6],'filled'); +plot([-1 1],[1/length(includedModels) 1/length(includedModels)],':','color','k','linewidth',0.5); +ylim([0 0.2]); ylabel('fraction of cue cells'); - -chance=1/90 * sum(bestModel(cueK&~behavior&~isnan(bestModel))<91)/sum(cueK&~behavior&~isnan(bestModel)); -plot([0 91],[chance chance],':','color','k','linewidth',0.5); -xlim([0 93]); -ylim([0 0.3]); -xticks(topModels(1:end-1)); +xlabel('correlation with value model'); yticks(0:0.1:0.3); -%xtickangle(45); -xticklabels(vnames(outliers(1:end-1))); -xtickangle(45); category=9*ones(length(improvement),1); -category(cueK&~behavior)=8; %cue neurons -%categorize neurons by their best model if it's in the top models -for bm=1:length(topModels) -category(cueK&~behavior&bestModel==topModels(bm))=bm; %cue, odor -end +category(cueK&~behavior)=4; %cue neurons, other odor coding schemes +category(cueK&~behavior&bestModel==1)=1; %value +category(cueK&~behavior&ismember(bestModel,valRank(2:17)))=2; %value like +category(cueK&~behavior&bestModel==max(bestModel))=3; %untuned + +% visual depiction of shuffles +colormap('summer'); +subplot(5,1,4); +imagesc(1:length(includedModels),1:6,valueShuffles(includedModels,:)',[0 1]); +xticks([]); +yticks(1:6); +yticklabels({'CS+','CS+','CS50','CS50','CS-','CS-'}); +title('Value assigned to each odor in each shuffle'); +colorbar; +% visual depiction of shuffles, reordered by correlation +subplot(5,1,5); +includedShuffles = valueShuffles(includedModels,:)'; -%% activity of cue cell types and correlations +imagesc(1:length(includedModels),1:6,includedShuffles(:,valRank),[0 1]); +xticks([]); +yticks(1:6); +yticklabels({'CS+','CS+','CS50','CS50','CS-','CS-'}); +title('Value assigned to each odor in each shuffle, ordered by similarity to value'); +colorbar; + +%% activity of cue cell types figure; map=redblue(256); -CL=[0 1]; + cueWindow=[0 2.5]; cueBins=binTimes>=cueWindow(1) & binTimes<=cueWindow(2); -cueCorrMatrix=NaN(12,12,totalNeurons); -cueCorrSets=NaN(6,6,totalNeurons); -setCorr=NaN(totalNeurons,3); - -for NN=1:totalNeurons - cueVector=NaN(sum(cueBins),12); - cue3Vector=NaN(sum(cueBins)*3,6); - v=0; - u=0; - for os=1:2 - for condition=1:6 - v=v+1; - cueVector(:,v)=dffPSTH6{condition,os}(NN,cueBins); - end - - end - %cueVector=cueVector-mean(cueVector,2); - corrmat=corrcoef(cueVector); - - cueCorrMatrix(:,:,NN)=corrmat; -end - - %heatmaps plotWindow=[-0.5 2.5]; plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); -sortWindow=[0 1.5]; +sortWindow=[0 2.5]; sortBins=binTimes>=sortWindow(1) & binTimes<=sortWindow(2); - - colors{1,1}=[0.1 0.6 0.2]; colors{2,1}=[0.4 0.1 0.4]; colors{3,1}=[0.3 0.3 0.3]; - -for ct=1 - sel=category==ct; -thisActivity=dffPSTH3{1,1}(sel,:); -cueResp=mean(thisActivity(:,sortBins),2); -[~,sortOrder]=sort(cueResp); pn=0; activity=[]; +thisActivity=dffPSTH3{1,1}; +cueResp=mean(thisActivity(:,sortBins),2); +cats={1,2,3}; + +regionCodeOpp=ones(totalNeurons,1); + +for ct=1:length(cats) +sel=ismember(category,cats{ct}); +sortcrit=[regionCodeOpp(sel) cueResp(sel)]; +[~,sortOrder]=sortrows(sortcrit,[1 2]); +pn=0; for cue=1:3 for os=1:2 pn=pn+1; - ax(1)=subplot(2,20,20+(ct-1)*6+pn); + ax(1)=subplot(3,15,(ct-1)*15+pn); colormap(ax(1),map); hold on; @@ -1012,11 +987,11 @@ - imagesc(binTimes(plotBins),[1 size(activity,1)],activity(:,plotBins),[-5 5]);%[min(min(cspActivity)) max(max(cspActivity))]); + imagesc(binTimes(plotBins),[1 size(activity,1)],activity(:,plotBins),[-3 3]);%[min(min(cspActivity)) max(max(cspActivity))]); ylim([0.5 size(activity,1)+0.5]) plot([0 0],[0.5 sum(sel)+0.5],'color',colors{cue},'linewidth',0.75); set(gca,'ytick',[]); - %if pn==1 ylabel('selective odor'); end + if pn==1 ylabel(sprintf('n=%d',sum(sel))); end if pn==1 xlabel('seconds from cue'); end if pn>1 xticks([]); end @@ -1026,106 +1001,61 @@ -end end +colors{1,1}=[0.1 0.6 0.2]; +colors{2,1}=[0.4 0.1 0.4]; +colors{3,1}=[0.3 0.3 0.3]; +directions=[1 -1]; +specs={'-','--'}; +diractivity = mean(cat(3,dffPSTH6{(1-1)*2+2,1}(:,plotBins),dffPSTH6{(1-1)*2+1,2}(:,plotBins)),3); +direction=sign(max(diractivity,[],2)-abs(min(diractivity,[],2))); - -pn=0; -activity=[]; -sel=ismember(category,2:6); -thisActivity=dffPSTH3{1,1}(sel,:); -cueResp=mean(thisActivity(:,sortBins),2); - -regionCodeOpp=8-category; -sortcrit=[regionCodeOpp(sel) cueResp]; -[~,sortOrder]=sortrows(sortcrit,[1 2]); -for cue=1:3 - for os=1:2 - - pn=pn+1; + for d=1:2 + subplot(6,3,3+(ct-1)*6+(d-1)*3); - ax(1)=subplot(2,20,30+(ct-1)*6+pn); - colormap(ax(1),map); hold on; + for cue=1:3 + for os=1:2 + + activity = dffPSTH6{(cue-1)*2+os,os}; + + %get values + psth=nanmean(activity(sel&direction==directions(d),plotBins)); + sem=nanste(activity(sel&direction==directions(d),plotBins),1); %calculate standard error of the mean + up=psth+sem; + down=psth-sem; + + %plotting + plot(binTimes(plotBins),psth,specs{os},'Color',colors{cue,1},'linewidth',0.75); + %patch([xvals,xvals(end:-1:1)],[up,down(end:-1:1)],regcolors{reg,1},'EdgeColor','none');alpha(0.2); + + + end + + end + if d==1 text(0.1,1.5,num2str(sum(sel&direction==directions(d)))); end + if d==2 text(0.1,-0.45,num2str(sum(sel&direction==directions(d)))); end + if d==2 + xlabel('seconds from odor onset'); + ylabel('z-score'); + else + xticks([]); + end - activity = dffPSTH3{cue,os}(sel,:); - activity = activity(sortOrder,:); - - - - imagesc(binTimes(plotBins),[1 size(activity,1)],activity(:,plotBins),[-5 5]);%[min(min(cspActivity)) max(max(cspActivity))]); - ylim([0.5 size(activity,1)+0.5]) - plot([0 0],[0.5 sum(sel)+0.5],'color',colors{cue},'linewidth',0.75); - set(gca,'ytick',[]); - %if pn==1 ylabel('selective odor'); end - if pn==1 xlabel('seconds from cue'); end - if pn>1 xticks([]); end - - - + %if reg>1 yticks([]); end + plot([0 0],[-1 5],'color','k','linewidth',0.75); + plot(plotWindow,[0 0],':','color','k','linewidth',0.5); + if d==1 & ct==1 axis([plotWindow -0.25 5]); end + if d==1 & ct==2 axis([plotWindow -0.25 5]); end + if d==2 axis([plotWindow -0.75 0.5]); end end - - - -end -titles={vnames{outliers(1:end-1)},'untuned','odor'}; -for ct=1:7 -ax(2)=subplot(4,7,7+ct); -sel=category==ct; -imagesc(nanmean(cueCorrMatrix(:,:,sel),3),CL); -%title(titles(ct)); -colormap(ax(2),'bone'); -xticks([]); -yticks([]); -end -plotWindow=[-0.5 2.5]; -plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); -cn=0; -for ct=1:7 -sel=category==ct; -if sum(sel)>0 - cn=cn+1; -activity=[]; -for cue=1:3 - for os=1:2 - activity = [activity dffPSTH3{cue,os}(sel,plotBins)]; - end end -activity=activity./max(abs(activity),[],2); -[coeff,score,~,~,explained]=pca(activity'); -specs={'-','--'}; -for comp=1 - cond=0; - if explained(comp)>5 - subplot(4,7,ct) - - hold on - for cue=1:3 - for os=1:2 - cond=cond+1; - plot(binTimes(plotBins),score((cond-1)*sum(plotBins)+1:cond*sum(plotBins),comp),specs{os},'color',colors{cue},'linewidth',1); - end - end - ylabel(sprintf('#%g (%g%%)',comp,round(explained(comp),1))); - yticks([]); - xticks([0 2]); - if comp==1 title(titles{cn}); end - -% plot([0 0],[min(score(:,comp)) max(score(:,comp))],':','color','k','linewidth',0.5); -% plot(plotWindow,[0 0],':','color','k','linewidth',0.5); - - end - xlim(cueWindow); -end -end -end - %% remove doubles (from merging) and edge of FOV rois function iscell = processROIs(iscell,stat,F) ROIpixels={};