diff --git a/Code/ephysMain.m b/Code/ephysMain.m index bb1347f..932399d 100644 --- a/Code/ephysMain.m +++ b/Code/ephysMain.m @@ -200,6 +200,21 @@ regleg{reg}=sprintf('%s (%d)',regions{reg},regionNeurons(reg)); end legend(regleg); +%% +figure; +subjects=unique(neuronSubject); +subjregn=zeros(length(subjects),length(regions)); +for subj=1:length(subjects) + for reg=1:length(regions) + subjregn(subj,reg)=sum(ismember(neuronSubject,subjects(subj))&ismember(neuronRegionOlf,regions(reg))); + end +end +heatmap(subjregn); +% colormap('summer'); +% yticks(1:length(subjects)); +% yticklabels(subjects); +% xticks(1:length(regions)); +% xticklabels(regions); %% analyze behavior subjects=unique(neuronSubject); @@ -468,6 +483,25 @@ text(5,0,sprintf('p = %g',round(p(1),2,'significant'))); end + +subplot(4,3,6) +hold on +for subject=1:length(subjects) + scatter(mean(anticipatoryLicks{2}(ismember(sessionSubject(includedSessions),subjects(subject)),:),2),mean(anticipatoryLicks{1}(ismember(sessionSubject(includedSessions),subjects(subject)),:),2),18,'filled'); +end + +plot([-2 10],[-2 10],':','color','k','linewidth',1); +axis([-2 10 -2 10]) +xticks([0 5 10]); +yticks([0 5 10]); +text(5,1,sprintf('combined')); + +data=[mean(anticipatoryLicks{2}(:,:),2);mean(anticipatoryLicks{1}(:,:),2)]; +gl=[zeros(length(data)/2,1);ones(length(data)/2,1)]; +sl=[sessionSubject(includedSessions);sessionSubject(includedSessions)]; +[p,tbl,stats]=anovan(data(:),{gl(:),sl(:)},'varnames',{'cue','subject'},'model','interaction','display','off'); +text(5,0,sprintf('p = %g',round(p(1),2,'significant'))); + %% plot example session trials figure; subplot(5,1,1); @@ -932,6 +966,7 @@ cueResp=mean(cspActivity(sel,cueBins),2); regions={'ALM','MOs','ACA','FRP','PL','ILA','ORB','CP','ACB','DP','TTd','AON','OLF'}; + regcolors{1,1}=[1 0.8 0]; regcolors{2,1}=[0.8 0.8 0.6]; regcolors{3,1}=[0.9 0.5 0.2]; @@ -996,6 +1031,175 @@ +end + +%% plot heatmaps of 3 main trial types PL only +map=redblue(256); +figure; +colormap(map); +neuronRegionOlf=neuronRegionAdj; +neuronRegionOlf(ismember(neuronRegionAdj,{'EPd','PIR'}))={'OLF'}; +sel=ismember(neuronRegionOlf,'PL'); + +%heatmaps +cspActivity=dffPSTH3{1,1}; %get the firing rates of neurons of interest +cueWindow=[0 1.5]; +cueBins=binTimes>=cueWindow(1) & binTimes<=cueWindow(2); +cueResp=mean(cspActivity(sel,cueBins),2); +% [~,cueResp]=max(abs(cspActivity(sel,cueBins)),[],2); +% cueResp=-cueResp; + +sortcrit=[ones(sum(sel),1) cueResp]; +[~,sortOrder]=sortrows(sortcrit,[1 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]; + +plotWindow=[-0.5 6]; +plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); +titles={'CS+ ','CS50 ','CS- '}; +sets={'set1','set2'}; +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; + subplot(1,6,pn); + hold on; + + imagesc(binTimes(plotBins),[1 sum(sel)],activity(:,plotBins),[-3 3]);%[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'); + end + end + + + +end + +%% PCA and correlations for PL only + +figure; + +plotWindow=[-0.5 6]; +plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); +sel=ismember(neuronRegionOlf,'PL'); +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:5 + cond=0; + if explained(comp)>5 + subplot(2,2,comp) + + 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([]); + if comp==1 xticks([0 2 4 6]); end + if comp>1 xticks([]); end + +% plot([0 0],[min(score(:,comp)) max(score(:,comp))],':','color','k','linewidth',0.5); +% plot(plotWindow,[0 0],':','color','k','linewidth',0.5); + xlim(plotWindow); + + patch([0 0 stimDuration stimDuration],[-1 1 1 -1],[0.6 0.3 0],'edgecolor','none');alpha(0.3); + plot([rewOnset rewOnset],[-1 1],'color',[0 0.6 0.3]); + end + +end + +%% plot mean activity for each region +figure; + +plotWindow=[-0.5 2.5]; +plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); +xvals=binTimes(plotBins); + +neuronRegionOlf=neuronRegionAdj; +neuronRegionOlf(ismember(neuronRegionAdj,{'EPd','PIR'}))={'OLF'}; +division=neuronRegionOlf; + +regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; +regcolors{1,1}=[1 0.8 0]; +regcolors{2,1}=[0.6 0.2 0.6]; +regcolors{3,1}=[1 0.6 1]; +regcolors{4,1}=[0.3 0.7 0.5]; + + +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},dffPSTH6{(1-1)*2+1,2}),3); +direction=sign(max(diractivity,[],2)-abs(min(diractivity,[],2))); +for reg=1:length(regions) + regsel=ismember(division,regions(reg)); + for d=1:2 + subplot(2,length(regions),reg+(d-1)*length(regions)); + + hold on; + for cue=1:3 + for os=1:2 + + activity = dffPSTH6{(cue-1)*2+os,os}; + + %get values + psth=nanmean(activity(regsel&direction==directions(d),plotBins)); + sem=nanste(activity(regsel&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(regsel&direction==directions(d)))); end + if d==2 text(0.1,-0.45,num2str(sum(regsel&direction==directions(d)))); end + if reg==1 & d==2 + xlabel('seconds from odor onset'); + ylabel('z-score'); + else + xticks([]); + end + + if reg>1 yticks([]); end + plot([0 0],[-1 2],'color',colors{condition},'linewidth',0.75); + plot(plotWindow,[0 0],':','color','k','linewidth',0.5); + if d==1 axis([plotWindow -0.25 2]); end + if d==2 axis([plotWindow -0.6 0.2]); end + end + + + end @@ -1087,6 +1291,7 @@ category(predictive & behavior & cueK,1)=2; category(predictive & behavior & rewardK,1)=4; category(predictive & behavior & cueK & rewardK,1)=7; +%category(sum(improvement(:,2:4)>0.05,2)>0)=8; %region neuronRegionOlf=neuronRegionAdj; @@ -1277,9 +1482,12 @@ barData=barData(:,regInd); b=bar(barData','stacked'); - + %b(3).FaceColor=[0.4 0 0.5]; %odor + %b(2).FaceColor=[0.7 0.7 0.7]; %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=[]; @@ -1351,7 +1559,12 @@ barData=barData(:,reggrps); b=bar(barData','stacked'); + %b(3).FaceColor=[0.4 0 0.5]; %odor + %b(2).FaceColor=[0.7 0.7 0.7]; %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=[]; @@ -1636,7 +1849,7 @@ direction=sign(max(diractivity,[],2)-abs(min(diractivity,[],2))); activity = mean(cat(3,dffPSTH6{(condition-1)*2+2,1},dffPSTH6{(condition-1)*2+1,2}),3); for d=1:2 - subplot(2,6,condition+3+(d-1)*6); + subplot(4,6,condition+3+(d-1)*6); hold on; for reg=1:length(cats) regsel=ismember(category,cats(reg)); @@ -1670,15 +1883,176 @@ if d==1 axis([plotWindow -0.25 2.25]);yticks(0:1:2); end if d==2 axis([plotWindow -1 0.2]);yticks(-1:0.5:0); end %if condition==1 legend(regions); end + + if condition==2 + + for reg=1:3 + subplot(4,6,reg+3+12+(d-1)*6); + hold on; + regsel=ismember(category,cats(reg)); + if cats(reg)==4 regsel=category>3; end + + for r=1:2 + + activity = mean(cat(3,dffPSTH4{r+1,1},dffPSTH4{r+1,2}),3); + %get values + psth=nanmean(activity(regsel&direction==directions(d),plotBins)); + sem=nanste(activity(regsel&direction==directions(d),plotBins),1); %calculate standard error of the mean + up=psth+sem; + down=psth-sem; + + %plotting + plot(binTimes(plotBins),psth,'Color',catcolors{reg,1}/r,'linewidth',0.75); + patch([xvals,xvals(end:-1:1)],[up,down(end:-1:1)],catcolors{reg,1}/r,'EdgeColor','none');alpha(0.2); + + + end + if reg==1 & d==2 + xlabel('seconds from odor onset'); + ylabel('z-score'); + else + xticks([]); + end + + if reg>1 yticks([]); end + plot([2.5 2.5],[-2 4],'color','k','linewidth',0.75); + patch([0 0 stimDuration stimDuration],[-2 4 4 -2],colors{condition},'edgecolor','none');alpha(0.3); + plot(plotWindow,[0 0],':','color','k','linewidth',0.5); + if d==1 axis([plotWindow -0.25 2]);yticks(0:1:2); end + if d==2 axis([plotWindow -1 0.2]);yticks(-1:0.5:0); end + end + + + end end end +%% single probe data +addpath(genpath(fullfile(githubDir, 'spikes'))) %to load .npy files +probe2plot = 28; %28 %3 +probesel=ismember(neuronProbe,probe2plot); + +plotWindow=[-0.5 6]; +plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); +msize=25; +regions={'MOs','ACA','PL','ILA','DP','TTd'}; +regcolors{1,1}=[0.8 0.8 0.6]; +regcolors{2,1}=[0.9 0.5 0.2]; +regcolors{3,1}=[0.6 0.2 0.6]; +regcolors{4,1}=[0.5 0.4 0.7]; +regcolors{5,1}=[0.3 0.7 0.5]; +regcolors{6,1}=[0.1 0.8 0.8]; + +colors={[0 0.4 0.9];... %cue + [0.6 0 0.6];... + [0.9 0.2 0];... %lick + [0.7 0.7 0.7]}; +neuronXYZmm=neuronXYZ/1000; + +subplot(1,8,1); +hold on; +for ct=1:length(regions) + sel=probesel & ismember(neuronRegionOlf,regions(ct)); + offset=(rand(sum(sel),1)*10-5)/1000; + p1=scatter(abs(neuronXYZmm(sel,1))+offset,neuronXYZmm(sel,3)+offset,msize,regcolors{ct},'filled'); + p1.MarkerFaceAlpha=0.5; +end +xlabel('ML'); +ylabel('DV'); +legend(regions,'location','se'); +title(sprintf('Insertion %d', probe)); +axis([200 1200 -4500 -1700]/1000); +NN=0; +NP=0; +clusterTimes={}; +clusterRegion={}; +clusterDepth=[]; +for session=1:sessionNumber + sessionFolder=fullfile(direc,sessionSubject{session},sessionDate{session}); + 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); + probeNeur=0; + NP=NP+1; + if NP==probe2plot + 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 + shankN=NaN(shanks,1); + 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'))); + channelLocationsID = readNPY(fullfile(shankFolder,'channels.CCFregionIDs.npy')); + channelXYZ = readNPY(fullfile(shankFolder,'channels.ML_AP_DV_fromBregma.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')); + shankN(shank)=length(incClusters); + + %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(regions,noNums(1:min([3 length(noNums)]))))==1 + NN=NN+1; + clusterTimes{NN,1} = double(spikeTimes(spikeClusters==cluster))/30000; + clusterRegion{NN,1} = noNums(1:min([3 length(noNums)])); + clusterDepth(NN,1)=channelXYZ(clusterChannel(clusterNames==cluster),3); + end + end + end + end + end + end +end + +regions={'PL'}; +regcolors{1,1}=[0.6 0.2 0.6]; + +subplot(1,5,3:5) +hold on +NN=0; +[ds,depthSorted]=sort(clusterDepth,'descend'); +for reg=1:length(regions) + for cluster=1:length(clusterTimes) + clusterSorted=depthSorted(cluster); + if strcmp(clusterRegion{clusterSorted},regions{reg}) + [xx,yy]=rasterize(clusterTimes{clusterSorted}); + plot(xx,yy+NN,'color',regcolors{reg},'linewidth',0.25); + NN=NN+1; + end + end +end +set(gca,'ydir','rev'); +xlim([850 1000]); +xticks([850 910]); +ylim([0 NN]); %% region comparisons of cue, lick, and both cells for striatum and pfc figure; regions={'ACA','FRP','PL','ILA','ORB','CP','ACB'}; @@ -1712,7 +2086,12 @@ barData=barData(:,regInd); b=bar(barData','stacked'); + %b(3).FaceColor=[0.4 0 0.5]; %odor + %b(2).FaceColor=[0.7 0.7 0.7]; %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=[]; @@ -1784,7 +2163,12 @@ barData=barData(:,reggrps); b=bar(barData','stacked'); + %b(3).FaceColor=[0.4 0 0.5]; %odor + %b(2).FaceColor=[0.7 0.7 0.7]; %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/ephysValue.m b/Code/ephysValue.m index 124b951..6d76acd 100644 --- a/Code/ephysValue.m +++ b/Code/ephysValue.m @@ -4,7 +4,6 @@ addpath(genpath(fullfile(githubDir, 'steinmetz-et-al-2019'))) %for reduced rank direc = 'D:\GitHub\ottenheimer-et-al-2022\Neuropixels'; %where neuropixels data exists addpath(direc); - %% find all sessions subjectFolders = dir(direc); @@ -1000,7 +999,7 @@ tic -%permute all possible combinations of odors for meaning shuffle +%permute all possible combinations of odors for value shuffle allodors=1:6; perms6=NaN(90,6); first2=nchoosek(allodors,2); @@ -1153,7 +1152,12 @@ 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}; @@ -1180,7 +1184,8 @@ end toc -%% value coding analysis 90 + +%% value coding analysis 90 with untuned improvCutoff=0.02; overallCutoff=0.02; improvement=varExp(:,1)-varExp; @@ -1199,7 +1204,7 @@ %value specific variance %whatever is accounted for by value kernel -cueVariance=totalCueVariance - valueImp(:,[4:end-1]); +cueVariance=totalCueVariance - valueImp(:,[4:end]); [~,bestModel]=max(cueVariance,[],2); bestModel(isnan(cueVariance(:,1)))=NaN; %region @@ -1217,44 +1222,71 @@ 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 subplot(3,1,1); hold on -frequency=histcounts(bestModel(cueK&~behavior&~isnan(bestModel)),0.5:1:90.5,'normalization','probability'); -outliers=frequency>0.03; +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),frequency,'facecolor',[0.6 0.6 0.6],'linewidth',0.01); +b=bar([1:max(bestModel)-1 max(bestModel)+1],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]; %meaning +b.CData(topModels(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 ylabel('fraction of cue cells'); -chance=1/90; +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 91]); -ylim([0 0.375]); -xticks(topModels); +xlim([0 93]); +ylim([0 0.35]); +xticks(topModels(1:end-1)); yticks(0:0.1:0.3); %xtickangle(45); %xticklabels(vnames(outliers)); -category=8*ones(length(improvement),1); -category(cueK&~behavior)=7; %cue neurons +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 -%% get proportions of value and meaning and compare between regions +% visual depiction of shuffles +colormap('summer'); +subplot(5,1,3); +imagesc(1:45,1:6,shuffVals(1:45,:)',[0 0.5]); +xticks(topModels); +yticks(1:6); +yticklabels({'CS+','CS+','CS50','CS50','CS-','CS-'}); +title('Value assigned to each odor in each shuffle'); + +subplot(5,1,4); +imagesc(46:90,1:6,shuffVals(46:end,:)',[0 0.5]); +xticks(topModels); +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); +colorbar; + +%% get proportions of value and trial type and compare between regions regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; reggrps=[1:3]; @@ -1268,16 +1300,16 @@ 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]}; -titles={'value','meaning (non-value)','non-meaning'}; -cats={1,2:6,7}; -for ct=1:3 +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 %region barData=zeros(2,5); - subplot(4,3,ct); + subplot(4,4,ct); hold on; sel=ismember(category,cats{ct}); - othersel=ismember(category,1:7)&~sel; + othersel=ismember(category,1:8)&~sel; regionnumber=zeros(totalNeurons,1); for reg=1:length(allReg) regsel=ismember(neuronRegionOlf,allReg(reg)); @@ -1291,12 +1323,8 @@ barData=barData(:,regInd); b=bar(barData','stacked'); - %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]; + b(2).FaceColor=[0.85 0.85 0.85]; + b(1).FaceColor=catcolors{ct}; upperCI=[]; lowerCI=[]; estmean=[]; @@ -1334,7 +1362,7 @@ greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end - s=subplot(4,3,3+ct); + s=subplot(4,4,4+ct); imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); @@ -1353,10 +1381,10 @@ %region group barData=zeros(2,3); - subplot(4,3,6+ct); + subplot(4,4,8+ct); hold on; sel=ismember(category,cats{ct}); - othersel=ismember(category,1:7)&~sel; + othersel=ismember(category,1:8)&~sel; regionnumber=zeros(totalNeurons,1); for reg=1:4 regsel=ismember(neuronGroup,reg); @@ -1370,12 +1398,8 @@ 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}; %value - %b(4).FaceColor=[0 0 0.2]; %odor - % b(5).FaceColor=[0.9 0.2 0]; - %b(3).FaceColor=[1 1 1]; + b(1).FaceColor=catcolors{ct}; upperCI=[]; lowerCI=[]; estmean=[]; @@ -1413,7 +1437,7 @@ greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end - s=subplot(4,3,9+ct); + s=subplot(4,4,12+ct); imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); @@ -1433,7 +1457,7 @@ end -%% get proportions of value and meaning out of cue cells +%% get proportions of value and trial type out of cue cells regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON'}; reggrps=[1:3]; @@ -1447,19 +1471,19 @@ 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]}; -titles={'value','meaning (non-value)','non-meaning'}; -cats={1,2:6,7}; -for ct=1:3 +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 %region barData=zeros(2,5); - subplot(4,3,ct); + subplot(4,4,ct); hold on; sel=ismember(category,cats{ct}); - othersel=ismember(category,1:7)&~sel; + othersel=ismember(category,1:8)&~sel; regionnumber=zeros(totalNeurons,1); for reg=1:length(allReg) - regsel=ismember(neuronRegionOlf,allReg(reg))&category<8; + regsel=ismember(neuronRegionOlf,allReg(reg))&category<9; regionnumber(regsel)=reg; barData(1,reg)=sum(sel®sel)/sum(regsel); barData(2,reg)=sum(othersel®sel)/sum(regsel); @@ -1470,13 +1494,13 @@ barData=barData(:,regInd); b=bar(barData','stacked'); - b(2).FaceColor=[0.85 0.85 0.85]; %meaning - b(1).FaceColor=catcolors{ct}; %value + b(2).FaceColor=[0.85 0.85 0.85]; + b(1).FaceColor=catcolors{ct}; upperCI=[]; lowerCI=[]; estmean=[]; for reg=1:length(regions) - regsel=find(ismember(neuronRegionOlf(category<8),regions(reg)),1); + regsel=find(ismember(neuronRegionOlf(category<9),regions(reg)),1); estmean(1,reg)=ypred(regsel); upperCI(1,reg)=ypredCI(regsel,2); lowerCI(1,reg)=ypredCI(regsel,1); @@ -1509,7 +1533,7 @@ greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end - s=subplot(4,3,3+ct); + s=subplot(4,4,4+ct); imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); @@ -1528,13 +1552,13 @@ %region group barData=zeros(2,3); - subplot(4,3,6+ct); + subplot(4,4,8+ct); hold on; sel=ismember(category,cats{ct}); - othersel=ismember(category,1:7)&~sel; + othersel=ismember(category,1:8)&~sel; regionnumber=zeros(totalNeurons,1); for reg=1:4 - regsel=ismember(neuronGroup,reg)&category<8; + regsel=ismember(neuronGroup,reg)&category<9; regionnumber(regsel)=reg; barData(1,reg)=sum(sel®sel)/sum(regsel); barData(2,reg)=sum(othersel®sel)/sum(regsel); @@ -1547,12 +1571,12 @@ barData=barData(:,reggrps); b=bar(barData','stacked'); b(2).FaceColor=[0.85 0.85 0.85]; - b(1).FaceColor=catcolors{ct}; %value + b(1).FaceColor=catcolors{ct}; upperCI=[]; lowerCI=[]; estmean=[]; for reg=1:length(reggrps) - regsel=find(ismember(neuronGroup(category<8),reggrps(reg)),1); + regsel=find(ismember(neuronGroup(category<9),reggrps(reg)),1); estmean(1,reg)=ypred(regsel); upperCI(1,reg)=ypredCI(regsel,2); lowerCI(1,reg)=ypredCI(regsel,1); @@ -1585,7 +1609,7 @@ greater(r1,r2)=lowerCI(1,r1)>upperCI(1,r2); end end - s=subplot(4,3,9+ct); + s=subplot(4,4,12+ct); imagesc(frac .* greater, [0 ceil(max(frac,[],'all')*20)/20]); colormap(s,customMap); @@ -1607,7 +1631,7 @@ %% value models schematic examplePerms=[1 33]; -exampleNeurons=[9 11 14 22 25 30]; +exampleNeurons=[9]; plotWindow=[0 2.5]; plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); @@ -1636,7 +1660,7 @@ for n=1:length(exampleNeurons) nn=entries(exampleNeurons(n)); ta=[]; - subplot(length(exampleNeurons),3,(n-1)*3+1); + subplot(length(exampleNeurons),5,(n-1)*5+1); %subplot(6,6,n); for os=1:2 @@ -1753,6 +1777,7 @@ for cue=1:6 originalValues=trialValues{session}(originalOrder==cue); newValue(newOrder==cue)=mean(originalValues(randsample(sum(originalOrder==cue),sum(newOrder==cue),'true'))); + newCueValue(perms6all(perm,cue))=mean(originalValues(randsample(sum(originalOrder==cue),sum(newOrder==cue),'true'))); end discreteUnfixed.values={newValue,ones(length(lickTimes),1),ones(length(boutStart),1)}; @@ -1770,11 +1795,12 @@ kernel=bR3(:,1:components)*fitK; pred(test)=A3(test,:)*kernel; end - ve = 1- var(y-pred)/var(y); + ve(n,permn) = 1- var(y-pred)/var(y); trialPrediction=reshape(pred,[binsPerTrial length(pred)/binsPerTrial])'; %3 conditions + examplePSTH={}; selections={}; selections{1,1}=csp1; %cs+ selections{2,1}=csf1; %cs50 @@ -1792,15 +1818,19 @@ end - subplot(length(exampleNeurons),3,(n-1)*3+1+permn); + subplot(length(exampleNeurons),5,(n-1)*5+1+permn); + hold on; + o=0; + for cue=1:3 for os=1:2 + o=o+1; - hold on; - for cue=1:3 - activity=examplePSTH{cue,os}(plotBinsGLM); - plot(binTimesGLMTrl(plotBinsGLM),activity,specs{os},'linewidth',1,'color',colors{cue}); - + + +% activity=examplePSTH{cue,os}(plotBinsGLM); +% plot(binTimesGLMTrl(plotBinsGLM),activity,specs{os},'linewidth',1,'color',colors{cue}); + plot(binTimesGLMTrl(plotBinsGLM),kernel(1:sum(plotBinsGLM))*newCueValue(o),specs{os},'linewidth',1,'color',colors{cue}); end end @@ -1889,7 +1919,6 @@ 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 @@ -1918,7 +1947,7 @@ pn=pn+1; - ax(1)=subplot(2,20,30+(ct-1)*6+pn); + ax(1)=subplot(2,20,27+(ct-1)*6+pn); colormap(ax(1),map); hold on; @@ -1932,7 +1961,6 @@ 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 @@ -1944,11 +1972,50 @@ end +pn=0; +activity=[]; +sel=ismember(category,7:8); +thisActivity=PSTH3{1,1}(sel,:); +cueResp=mean(thisActivity(:,sortBins),2); -titles={vnames{outliers},'non-meaning'}; +regionCodeOpp=8-category; +sortcrit=[regionCodeOpp(sel) cueResp]; +[~,sortOrder]=sortrows(sortcrit,[1 2]); +for cue=1:3 + for os=1:2 + + pn=pn+1; + + ax(1)=subplot(2,20,34+(ct-1)*6+pn); + colormap(ax(1),map); + hold on; + + + activity = PSTH3{cue,os}(sel,:); + activity = activity(sortOrder,:); + -for ct=1:7 -ax(2)=subplot(4,7,7+ct); + + 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 + + +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)); @@ -1959,7 +2026,7 @@ plotWindow=[-0.5 2.5]; plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); -for ct=1:7 +for ct=1:8 sel=category==ct; activity=[]; for cue=1:3 @@ -1973,7 +2040,7 @@ for comp=1 cond=0; if explained(comp)>5 - subplot(4,7,ct) + subplot(4,8,ct) hold on for cue=1:3 @@ -1995,8 +2062,321 @@ end end +%% plot PC1 value from each region +figure; + +cueWindow=[0 2.5]; +cueBins=binTimes>=cueWindow(1) & binTimes<=cueWindow(2); +neuronRegionOlf=neuronRegionAdj; +neuronRegionOlf(ismember(neuronRegionAdj,{'EPd','PIR'}))={'OLF'}; +division=neuronRegionOlf; +regions={'ALM','ACA','FRP','PL','ILA','ORB','DP','TTd','AON','CP','ACB'}; + +plotWindow=[-0.5 2.5]; +plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); +for reg=1:length(regions) + sel=ismember(division,regions(reg)) & category==1; +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={'-','--'}; +comp=1; +if strcmp(regions(reg),'ACA') comp=2; end +cond=0; +if explained(comp)>5 + subplot(4,9,reg) + + 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]); + title(regions{reg}); + + % 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); +xlabel(sprintf('n = %d',sum(sel))); + + + +end + +division=neuronGroup; +for reg=1:4 + sel=ismember(division,reg) & category==1; +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,9,18+reg) + + 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]); + xlabel(sprintf('n = %d',sum(sel))); + if comp==1 title(regionGroupNames{reg}); 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 + +%% 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; @@ -2026,6 +2406,238 @@ end +%% coding dimension for CS+/CS50/CS- +figure; +%categories for subselection +sels={}; +sels{1}=ismember(category,8); +sels{2}=ismember(category,2:6); +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]}; + +numstraps=5000; +amin=-0.1; +amax=1.2; +cueWindow=[0 2.5]; +cueBins=binTimes>=cueWindow(1) & binTimes<=cueWindow(2); + +plotWindow=[-1 10]; +plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); +plotBinTimes=binTimes(plotBins); +baseBins=plotBinTimes<0; + +%normalize activity from 0-1 +dffPSTH6n={}; +dffPSTH3n={}; +for session=1 + dffPSTH6all=[]; + for condition=1:6 + dffPSTH6all=cat(2,dffPSTH6all,PSTH6{condition,session}(:,cueBins)); + end +end +for session=1:2 + for condition=1:6 + dffPSTH6n{condition,session}=PSTH6{condition,session}./max(abs(dffPSTH6all),[],2); + end +end + +%get activity for dimensions +csmActivity=dffPSTH6n{5,1}(:,cueBins); +cspActivity=dffPSTH6n{1,1}(:,cueBins); +csfActivity=dffPSTH6n{3,1}(:,cueBins); + +bigBinSize=5; +bigBins=floor(sum(cueBins)/bigBinSize); +bsf=0; +for bb=1:bigBins + csmActivityB(:,bb)=mean(csmActivity(:,bsf+1:bsf+bigBinSize),2); + cspActivityB(:,bb)=mean(cspActivity(:,bsf+1:bsf+bigBinSize),2); + csfActivityB(:,bb)=mean(csfActivity(:,bsf+1:bsf+bigBinSize),2); + + bsf=bsf+bigBinSize; +end + +[~,PMind]=max(abs(cspActivityB-csmActivityB),[],2,'linear'); +[~,PFind]=max(abs(cspActivityB-csfActivityB),[],2,'linear'); +[~,FMind]=max(abs(csfActivityB-csmActivityB),[],2,'linear'); + +differenceVectorPM=cspActivityB(PMind)-csmActivityB(PMind); +differenceVectorPF=cspActivityB(PFind)-csfActivityB(PFind); +differenceVectorFM=csfActivityB(FMind)-csmActivityB(FMind); + + +sproj={}; +oproj={}; +disttime={}; + +%for getting 0 and 1, use z-score +csmActivity=PSTH6{5,1}(:,cueBins); +cspActivity=PSTH6{1,1}(:,cueBins); +csfActivity=PSTH6{3,1}(:,cueBins); + +bsf=0; +for bb=1:bigBins + csmActivityB(:,bb)=mean(csmActivity(:,bsf+1:bsf+bigBinSize),2); + cspActivityB(:,bb)=mean(cspActivity(:,bsf+1:bsf+bigBinSize),2); + csfActivityB(:,bb)=mean(csfActivity(:,bsf+1:bsf+bigBinSize),2); + + bsf=bsf+bigBinSize; +end + +projCorrStrp=NaN(numstraps,length(sels)); +baseCSMdist=NaN(numstraps,length(sels)); +cspfangle=NaN(numstraps,length(sels)); +for cl=1:length(sels) + sel=sels{cl}; + + PF1=differenceVectorPF(sel)' * cspActivityB(PFind(sel)); + PF0=differenceVectorPF(sel)' * csfActivityB(PFind(sel)); + PM1=differenceVectorPM(sel)' * cspActivityB(PMind(sel)); + PM0=differenceVectorPM(sel)' * csmActivityB(PMind(sel)); + FM1=differenceVectorFM(sel)' * csfActivityB(FMind(sel)); + FM0=differenceVectorFM(sel)' * csmActivityB(FMind(sel)); + + for cue=1:3 + sproj{cue,1}=(differenceVectorPM(sel)' * PSTH6{(cue-1)*2+2,1}(sel,:) - PM0) ./ (PM1 - PM0); + sproj{cue,2}=(differenceVectorPF(sel)' * PSTH6{(cue-1)*2+2,1}(sel,:) - PF0) ./ (PF1 - PF0); + sproj{cue,3}=(differenceVectorFM(sel)' * PSTH6{(cue-1)*2+2,1}(sel,:) - FM0) ./ (FM1 - FM0); + + oproj{cue,1}=(differenceVectorPM(sel)' * PSTH3{cue,2}(sel,:) - PM0) ./ (PM1 - PM0); + oproj{cue,2}=(differenceVectorPF(sel)' * PSTH3{cue,2}(sel,:) - PF0) ./ (PF1 - PF0); + oproj{cue,3}=(differenceVectorFM(sel)' * PSTH3{cue,2}(sel,:) - FM0) ./ (FM1 - FM0); + end + + subplot(3,4,1+(cl-1)*4); + hold on; + for cue=1:3 + plot(sproj{cue,1}(plotBins),sproj{cue,3}(plotBins),':','linewidth',1,'color',colors{cue}); + plot(sproj{cue,1}(cueBins),sproj{cue,3}(cueBins),'linewidth',1,'color',colors{cue}); + end + axis([amin amax amin amax amin amax]); + if cl==1 + title('same odor set'); + xlabel('CS- ---> CS+'); + ylabel('CS- ---> CS50'); + end + xticks([0 1]); + yticks([0 1]); + + 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}(cueBins),oproj{cue,3}(cueBins),'linewidth',1,'color',colors{cue}); + + end + axis([amin amax amin amax amin amax]); + xticks([]); + yticks([]); + if cl==1 title('other odor set'); end + + + %bootstrap a statistic + incNeur=find(sel); + + for strap=1:numstraps + ns=incNeur(randsample(length(incNeur),length(incNeur),'true')); + PF1=differenceVectorPF(ns)' * cspActivityB(PFind(ns)); + PF0=differenceVectorPF(ns)' * csfActivityB(PFind(ns)); + PM1=differenceVectorPM(ns)' * cspActivityB(PMind(ns)); + PM0=differenceVectorPM(ns)' * csmActivityB(PMind(ns)); + FM1=differenceVectorFM(ns)' * csfActivityB(FMind(ns)); + FM0=differenceVectorFM(ns)' * csmActivityB(FMind(ns)); + + for cue=1:3 + sproj{cue,1}=(differenceVectorPM(ns)' * PSTH6{(cue-1)*2+2,1}(ns,cueBins) - PM0) ./ (PM1 - PM0); + sproj{cue,2}=(differenceVectorPF(ns)' * PSTH6{(cue-1)*2+2,1}(ns,cueBins) - PF0) ./ (PF1 - PF0); + sproj{cue,3}=(differenceVectorFM(ns)' * PSTH6{(cue-1)*2+2,1}(ns,cueBins) - FM0) ./ (FM1 - FM0); + + oproj{cue,1}=(differenceVectorPM(ns)' * PSTH3{cue,2}(ns,cueBins) - PM0) ./ (PM1 - PM0); + oproj{cue,2}=(differenceVectorPF(ns)' * PSTH3{cue,2}(ns,cueBins) - PF0) ./ (PF1 - PF0); + oproj{cue,3}=(differenceVectorFM(ns)' * PSTH3{cue,2}(ns,cueBins) - FM0) ./ (FM1 - FM0); + end + + %only use CS- subtracted axes for correlation + sprojall=cat(1,sproj{:,[1 3]}); + oprojall=cat(1,oproj{:,[1 3]}); + projCorrStrp(strap,cl)=corr(sprojall(:),oprojall(:)); + + %baseline distance from CS- + basex=mean([sproj{1,1}(baseBins) sproj{2,1}(baseBins) sproj{3,1}(baseBins)]); + basey=mean([sproj{1,3}(baseBins) sproj{2,3}(baseBins) sproj{3,3}(baseBins)]); + baseCSMdist(strap,cl)=sqrt(basex^2+basey^2); + + %angle between CS+ at peak along (CS-/CS+) and CS50 at peak along (CS-/CS50) + [cspx,bin]=max(sproj{1,1}); %CSP + cspy=sproj{1,3}(bin); + [csfy,bin]=max(sproj{2,3}); %CSF + csfx=sproj{2,1}(bin); + v_1 = [cspx,cspy,0] - [basex,basey,0]; + v_2 = [csfx,csfy,0] - [basex,basey,0]; + cspfangle(strap,cl) = atan2(norm(cross(v_1, v_2)), dot(v_1, v_2)) * 180/pi; %angle + + end + + + + +end + +subplot(2,4,7); +hold on; +for c=1:3 +histogram(baseCSMdist(:,c),0.1:0.01:1,'facecolor',vcolors{c},'edgecolor',vcolors{c},'orientation','horizontal','normalization','probability'); +end + +ylabel('dist from cs-'); +xlabel('frequency'); +legend('odor','trial type','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); +text(0.1,0.8,sprintf('p = %g',round(bootp,2,'significant'))); +ylim([0 0.8]); + +subplot(2,4,8); +hold on; +for c=1:3 +histogram(cspfangle(:,c),0:1:110,'facecolor',vcolors{c},'edgecolor',vcolors{c},'orientation','horizontal','normalization','probability'); +end + +ylabel('CS+, CS50 angle'); +xlabel('frequency'); +legend('odor','trial type','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]); + + +subplot(2,4,4); +hold on; +for c=1:3 +histogram(projCorrStrp(:,c),0.8:0.001:1,'facecolor',vcolors{c},'edgecolor',vcolors{c},'orientation','horizontal','normalization','probability'); +end + +ylabel('same/other correlation'); +xlabel('frequency'); +legend('odor','trial type','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); +text(0.15,.97,sprintf('p = %g',round(bootp,3,'significant'))); +bootp=(sum(sum(projCorrStrp(:,1)>=projCorrStrp(:,2)'))+1)/(numstraps^2+1); +text(0.2,.93,sprintf('p = %g',round(bootp,3,'significant'))); +ylim([0.9 1]); +yticks([0.9 0.95 1]); + %% perform regression for shuffling trial values for value cells alreadyRun=true; @@ -2147,7 +2759,7 @@ hold on; sel=category==1 & bestValue==1; othervalue=ismember(category,1)&~sel; - othersel=ismember(category,2:7); + othersel=ismember(category,2:8); regionnumber=zeros(totalNeurons,1); for reg=1:length(allReg) regsel=ismember(neuronRegionOlf,allReg(reg)); @@ -2214,7 +2826,7 @@ xtickangle(45); yticks(1:length(regions)); yticklabels(regions); - title(titles{ct}); + %title(titles{ct}); if ct==1 ylabel('increase in this region...'); xlabel('over this region'); @@ -2281,7 +2893,7 @@ xtickangle(45); yticks(1:length(reggrps)); yticklabels(regionGroupNames(reggrps)); - title(titles{ct}); + %title(titles{ct}); if ct==1 ylabel('increase in this region...'); xlabel('over this region'); @@ -2313,7 +2925,7 @@ hold on; sel=category==1 & bestValue==1; othervalue=ismember(category,1)&~sel; - othersel=ismember(category,2:7); + othersel=ismember(category,2:8); regionnumber=zeros(totalNeurons,1); for reg=1:length(allReg) regsel=ismember(neuronRegionOlf,allReg(reg))&category==1; @@ -2380,7 +2992,7 @@ xtickangle(45); yticks(1:length(regions)); yticklabels(regions); - title(titles{ct}); + %title(titles{ct}); if ct==1 ylabel('increase in this region...'); xlabel('over this region'); @@ -2447,14 +3059,13 @@ xtickangle(45); yticks(1:length(reggrps)); yticklabels(regionGroupNames(reggrps)); - title(titles{ct}); + %title(titles{ct}); if ct==1 ylabel('increase in this region...'); xlabel('over this region'); end end - %% value models schematic examplePerms=[1 1]; exampleNeurons=[5 20 50 72 111 141]; @@ -2674,9 +3285,9 @@ sels{1}=category==1 & bestValue==1; sels{2}=category==1 & bestValue==2; sels{3}=ismember(category,2:6); -sels{4}=ismember(category,7); +sels{4}=ismember(category,7:8); sels{5}=behavior; -cattitles={'value (history)','value (non-history)','meaning (non-value)','non-meaning','lick, cue+lick'}; +cattitles={'value (history)','value (non-history)','trial type','other cue','lick, cue+lick'}; valcolors{1,1}=[0.2 0.2 0.2]; valcolors{2,1}=[0.5 0.5 0.5]; @@ -2870,13 +3481,13 @@ end -%% get proportions of value and meaning and compare between regions for striatum +%% get proportions of value and trial type 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','meaning (non-value)','non-meaning'}; +titles={'value','trial type (non-value)','non-trial type'}; cats={1,2:6,7}; for ct=1:3 %region @@ -2884,7 +3495,7 @@ subplot(4,3,ct); hold on; sel=ismember(category,cats{ct}); - othersel=ismember(category,1:7)&~sel; + othersel=ismember(category,1:8)&~sel; regionnumber=zeros(totalNeurons,1); for reg=1:length(regions) regsel=ismember(neuronRegionOlf,regions(reg)); @@ -2897,12 +3508,8 @@ [ypred,ypredCI]=predict(lm,'conditional',false); b=bar(barData','stacked'); - %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]; + b(2).FaceColor=[0.85 0.85 0.85]; + b(1).FaceColor=catcolors{ct}; upperCI=[]; lowerCI=[]; estmean=[]; @@ -2962,7 +3569,7 @@ subplot(4,3,6+ct); hold on; sel=ismember(category,cats{ct}); - othersel=ismember(category,1:7)&~sel; + othersel=ismember(category,1:8)&~sel; regionnumber=zeros(totalNeurons,1); regs=[2 4]; for reg=1:2 @@ -2976,8 +3583,12 @@ [ypred,ypredCI]=predict(lm,'conditional',false); 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}; %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/imagingOdorSets.m b/Code/imagingOdorSets.m index 7fb4342..8afab50 100644 --- a/Code/imagingOdorSets.m +++ b/Code/imagingOdorSets.m @@ -155,6 +155,7 @@ end figure; +subplot(1,4,1); for session=1:length(sessions) hold on; errorbar([1 2],mean(anticipatoryLicks{1,1}),nanste(anticipatoryLicks{1,1},1),'linewidth',1.5,'color',[0.1 0.6 0.2]); @@ -171,7 +172,21 @@ xlim([0.75 2.25]); end +subplot(1,2,2) +hold on +for subject=1:length(anticipatoryLicks{1}) + scatter(mean(anticipatoryLicks{2}(subject,:)),mean(anticipatoryLicks{1}(subject,:)),18,'filled'); +end + +plot([-2 10],[-2 10],':','color','k','linewidth',1); +axis([-2 10 -2 10]) +xticks([0 5 10]); +yticks([0 5 10]); +text(5,1,sprintf('combined')); +data=[mean(anticipatoryLicks{2}(:,:),2);mean(anticipatoryLicks{1}(:,:),2)]; +[~,p,~,stats]=ttest(mean(anticipatoryLicks{2}(:,:),2),mean(anticipatoryLicks{1}(:,:),2)); +text(5,0,sprintf('p = %g',round(p(1),2,'significant'))); %% count neurons totalNeurons=0; @@ -345,6 +360,7 @@ % medianF=medfilt1(double(correctedF),7); decSpks=spks(roiNum,firstFrame(mouse,session):lastFrame(mouse,session)); + %decSpks=smooth(decSpks,5)'; %to account for error in estimating deconvolved spike frame %smooth with the same filter used for ephys spikeRateSmooth=NaN(1,length(decSpks)); @@ -779,7 +795,265 @@ end +%% PCA +figure; + +plotWindow=[-0.5 6]; +plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); +sel=true(totalNeurons,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:5 + cond=0; + if explained(comp)>5 + subplot(2,2,comp) + + 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([]); + if comp==1 xticks([0 2 4 6]); end + if comp>1 xticks([]); end + +% plot([0 0],[min(score(:,comp)) max(score(:,comp))],':','color','k','linewidth',0.5); +% plot(plotWindow,[0 0],':','color','k','linewidth',0.5); + xlim(plotWindow); + + patch([0 0 stimDuration stimDuration],[-1 1 1 -1],[0.6 0.3 0],'edgecolor','none');alpha(0.3); + plot([rewOnset rewOnset],[-1 1],'color',[0 0.6 0.3]); + end + +end + + + +%% coding dimension for CS+/CS50/CS- +figure; +%categories for subselection +sels={}; +sels{1}=category{1}==1; +sels{2}=category{1}==3; + +numstraps=5000; +amin=-0.25; +amax=1.75; +cueWindow=[0 2.5]; +cueBins=binTimes>=cueWindow(1) & binTimes<=cueWindow(2); + +plotWindow=[-1 10]; +plotBins=binTimes>=plotWindow(1) & binTimes<=plotWindow(2); +plotBinTimes=binTimes(plotBins); +baseBins=plotBinTimes<0; + +%normalize activity from 0-1 +dffPSTH6n={}; +dffPSTH3n={}; +for session=1:2 + dffPSTH6all=[]; + for condition=1:6 + dffPSTH6all=cat(2,dffPSTH6all,dffPSTH6{condition,session}); + end + for condition=1:6 + dffPSTH6n{condition,session}=dffPSTH6{condition,session}./max(abs(dffPSTH6all),[],2); + end + + dffPSTH3all=[]; + for condition=1:3 + dffPSTH3all=cat(2,dffPSTH3all,dffPSTH3{condition,session}); + end + for condition=1:3 + dffPSTH3n{condition,session}=dffPSTH3{condition,session}./max(abs(dffPSTH3all),[],2); + end +end + +%get activity for dimensions +csmActivity=dffPSTH6n{5,1}(:,cueBins); +cspActivity=dffPSTH6n{1,1}(:,cueBins); +csfActivity=dffPSTH6n{3,1}(:,cueBins); + +bigBinSize=5; +bigBins=floor(sum(cueBins)/bigBinSize); +bsf=0; +for bb=1:bigBins + csmActivityB(:,bb)=mean(csmActivity(:,bsf+1:bsf+bigBinSize),2); + cspActivityB(:,bb)=mean(cspActivity(:,bsf+1:bsf+bigBinSize),2); + csfActivityB(:,bb)=mean(csfActivity(:,bsf+1:bsf+bigBinSize),2); + + bsf=bsf+bigBinSize; +end + +%normalize activity from -1 to 1 +[~,PMind]=max(abs(cspActivityB-csmActivityB),[],2,'linear'); +[~,PFind]=max(abs(cspActivityB-csfActivityB),[],2,'linear'); +[~,FMind]=max(abs(csfActivityB-csmActivityB),[],2,'linear'); + +differenceVectorPM=cspActivityB(PMind)-csmActivityB(PMind); +differenceVectorPF=cspActivityB(PFind)-csfActivityB(PFind); +differenceVectorFM=csfActivityB(FMind)-csmActivityB(FMind); + + +sproj={}; +oproj={}; +disttime={}; + +%for getting 0 and 1, use z-score +csmActivity=dffPSTH6{5,1}(:,cueBins); +cspActivity=dffPSTH6{1,1}(:,cueBins); +csfActivity=dffPSTH6{3,1}(:,cueBins); + +bsf=0; +for bb=1:bigBins + csmActivityB(:,bb)=mean(csmActivity(:,bsf+1:bsf+bigBinSize),2); + cspActivityB(:,bb)=mean(cspActivity(:,bsf+1:bsf+bigBinSize),2); + csfActivityB(:,bb)=mean(csfActivity(:,bsf+1:bsf+bigBinSize),2); + + bsf=bsf+bigBinSize; +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]; + +projCorrStrp=NaN(numstraps,length(sels)); +baseCSMdist=NaN(numstraps,length(sels)); +cspfangle=NaN(numstraps,length(sels)); +for cl=1:length(sels) + sel=sels{cl}; + + PF1=differenceVectorPF(sel)' * cspActivityB(PFind(sel)); + PF0=differenceVectorPF(sel)' * csfActivityB(PFind(sel)); + PM1=differenceVectorPM(sel)' * cspActivityB(PMind(sel)); + PM0=differenceVectorPM(sel)' * csmActivityB(PMind(sel)); + FM1=differenceVectorFM(sel)' * csfActivityB(FMind(sel)); + FM0=differenceVectorFM(sel)' * csmActivityB(FMind(sel)); + + for cue=1:3 + sproj{cue,1}=(differenceVectorPM(sel)' * dffPSTH6{(cue-1)*2+2,1}(sel,:) - PM0) ./ (PM1 - PM0); + sproj{cue,2}=(differenceVectorPF(sel)' * dffPSTH6{(cue-1)*2+2,1}(sel,:) - PF0) ./ (PF1 - PF0); + sproj{cue,3}=(differenceVectorFM(sel)' * dffPSTH6{(cue-1)*2+2,1}(sel,:) - FM0) ./ (FM1 - FM0); + + oproj{cue,1}=(differenceVectorPM(sel)' * dffPSTH3{cue,2}(sel,:) - PM0) ./ (PM1 - PM0); + oproj{cue,2}=(differenceVectorPF(sel)' * dffPSTH3{cue,2}(sel,:) - PF0) ./ (PF1 - PF0); + oproj{cue,3}=(differenceVectorFM(sel)' * dffPSTH3{cue,2}(sel,:) - FM0) ./ (FM1 - FM0); + end + + subplot(3,4,1+(cl-1)*4); + hold on; + for cue=1:3 + plot(sproj{cue,1}(plotBins),sproj{cue,3}(plotBins),':','linewidth',1,'color',colors{cue}); + plot(sproj{cue,1}(cueBins),sproj{cue,3}(cueBins),'linewidth',1,'color',colors{cue}); + + end + %view(45,0); + axis([amin amax amin amax amin amax]); + if cl==1 + title('same odor set'); + xlabel('CS- ---> CS+'); + ylabel('CS- ---> CS50'); + end + xticks([0 1]); + yticks([0 1]); + + 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}(cueBins),oproj{cue,3}(cueBins),'linewidth',1,'color',colors{cue}); + + end + axis([amin amax amin amax amin amax]); + xticks([]); + yticks([]); + if cl==1 title('other odor set'); end + + + %bootstrap a statistic + incNeur=find(sel); + + for strap=1:numstraps + ns=incNeur(randsample(length(incNeur),length(incNeur),'true')); + PF1=differenceVectorPF(ns)' * cspActivityB(PFind(ns)); + PF0=differenceVectorPF(ns)' * csfActivityB(PFind(ns)); + PM1=differenceVectorPM(ns)' * cspActivityB(PMind(ns)); + PM0=differenceVectorPM(ns)' * csmActivityB(PMind(ns)); + FM1=differenceVectorFM(ns)' * csfActivityB(FMind(ns)); + FM0=differenceVectorFM(ns)' * csmActivityB(FMind(ns)); + + for cue=1:3 + sproj{cue,1}=(differenceVectorPM(ns)' * dffPSTH6{(cue-1)*2+2,1}(ns,cueBins) - PM0) ./ (PM1 - PM0); + sproj{cue,2}=(differenceVectorPF(ns)' * dffPSTH6{(cue-1)*2+2,1}(ns,cueBins) - PF0) ./ (PF1 - PF0); + sproj{cue,3}=(differenceVectorFM(ns)' * dffPSTH6{(cue-1)*2+2,1}(ns,cueBins) - FM0) ./ (FM1 - FM0); + + oproj{cue,1}=(differenceVectorPM(ns)' * dffPSTH3{cue,2}(ns,cueBins) - PM0) ./ (PM1 - PM0); + oproj{cue,2}=(differenceVectorPF(ns)' * dffPSTH3{cue,2}(ns,cueBins) - PF0) ./ (PF1 - PF0); + oproj{cue,3}=(differenceVectorFM(ns)' * dffPSTH3{cue,2}(ns,cueBins) - FM0) ./ (FM1 - FM0); + end + + %only use CS- subtracted axes for correlation + sprojall=cat(1,sproj{:,[1 3]}); + oprojall=cat(1,oproj{:,[1 3]}); + projCorrStrp(strap,cl)=corr(sprojall(:),oprojall(:)); + + %baseline distance from CS- + basex=mean([sproj{1,1}(baseBins) sproj{2,1}(baseBins) sproj{3,1}(baseBins)]); + basey=mean([sproj{1,3}(baseBins) sproj{2,3}(baseBins) sproj{3,3}(baseBins)]); + baseCSMdist(strap,cl)=sqrt(basex^2+basey^2); + + %angle between CS+ at peak along (CS-/CS+) and CS50 at peak along (CS-/CS50) + [cspx,bin]=max(sproj{1,1}); %CSP + cspy=sproj{1,3}(bin); + [csfy,bin]=max(sproj{2,3}); %CSF + csfx=sproj{2,1}(bin); + v_1 = [cspx,cspy,0] - [basex,basey,0]; + v_2 = [csfx,csfy,0] - [basex,basey,0]; + cspfangle(strap,cl) = atan2(norm(cross(v_1, v_2)), dot(v_1, v_2)) * 180/pi; %angle + + end + + + + +end + +subplot(3,5,4); +hold on; +for c=1:2 +histogram(cspfangle(:,c),-10:1:90,'facecolor',catcolors{c},'edgecolor',catcolors{c},'orientation','horizontal','normalization','probability'); +end +ylabel('CS+, CS50 angle'); +xlabel('frequency'); +bootp=(sum(sum(cspfangle(:,1)<=cspfangle(:,2)'))+1)/(numstraps^2+1); +text(0.01,50,sprintf('p = %g',round(bootp,2,'significant'))); +ylim([-10 90]); +yticks([0 45 90]); +plot([0 0.1],[0 0],':','color','k','linewidth',0.5); +xticks([0 0.1]); + +subplot(3,5,5); +hold on; +for c=1:2 +histogram(projCorrStrp(:,c),0.5:0.01:1,'facecolor',catcolors{c},'edgecolor',catcolors{c},'orientation','horizontal','normalization','probability'); +end +ylabel('same/other correlation'); +bootp=(sum(sum(projCorrStrp(:,1)>=projCorrStrp(:,2)'))+1)/(numstraps^2+1); +text(0.1,0.7,sprintf('p = %g',round(bootp,2,'significant'))); +ylim([0.5 1]); +yticks([0.5 1]); +xticks([0 0.4]); +xlim([0 0.4]); %% remove doubles (from merging) and edge of FOV rois function iscell = processROIs(iscell,stat,F) diff --git a/Code/imagingOdorSetsValue.m b/Code/imagingOdorSetsValue.m index 8d72f2d..e6a6277 100644 --- a/Code/imagingOdorSetsValue.m +++ b/Code/imagingOdorSetsValue.m @@ -1,4 +1,4 @@ -% Script to analyze value coding across day 3 of each odor set +% Script to analyze olfactory conditioning data across day 3 of each odor set githubDir = 'D:\GitHub'; addpath(genpath(fullfile(githubDir, 'npy-matlab'))) addpath(genpath(fullfile(githubDir, 'steinmetz-et-al-2019'))) @@ -559,6 +559,7 @@ % medianF=medfilt1(double(correctedF),7); decSpks=spks(roiNum,firstFrame(mouse,session):lastFrame(mouse,session)); + %decSpks=smooth(decSpks,5)'; %to account for error in estimating deconvolved spike frame %smooth with the same filter used for ephys spikeRateSmooth=NaN(1,length(decSpks)); @@ -599,6 +600,16 @@ end +% %4 conditions +% selections{1,1}=cue1; %cs+ +% selections{2,1}=cue2r; %reward+ +% selections{3,1}=cue2u; %reward- +% selections{4,1}=cue3; %cs- +% for condition=1:length(selections) +% sel = ismember(round(cue),round(selections{condition})); +% dffPSTH4{condition,session}(NN,:)=nanmean(binneddFF(sel,:)); +% end + %3 conditions, odd and even trials selections={}; selections{1,1}=cue1; %cs+ @@ -644,6 +655,9 @@ includedActivitySh=binneddFF(shOrder,includedBins)'; activityVector=includedActivity(:); activityVectorSh=includedActivitySh(:); + %normalize from 0 to 1 +% activityVectorSh=activityVectorSh/std(activityVectorSh); +% activityVector=activityVector/std(activityVector); Yall{NS,1}(:,roi)=activityVector; YallSh{NS,1}(:,roi)=activityVectorSh; @@ -761,6 +775,9 @@ 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; @@ -795,7 +812,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 @@ -821,6 +843,11 @@ 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}; @@ -846,11 +873,12 @@ varExpValueSh(:,3)=varExp(:,5); -%% value coding analysis 90 + +%% value coding analysis 90 with untuned improvCutoff=0.02; overallCutoff=0.02; -improvement=-(varExp-varExp(:,1)); +improvement=varExp(:,1)-varExp; %get proportions predictive = varExp(:,1)>overallCutoff; %more than 2% of variance explained @@ -866,44 +894,52 @@ %value specific variance %whatever is accounted for by value kernel -cueVariance=totalCueVariance - valueImp(:,[4:end-1]); +cueVariance=totalCueVariance - valueImp(:,[4:end]); [~,bestModel]=max(cueVariance,[],2); +bestModel(isnan(cueVariance(:,1)))=NaN; 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 subplot(3,1,1); hold on -frequency=histcounts(bestModel(cueK&~behavior),0.5:1:90.5,'normalization','probability'); -outliers=frequency>0.04; +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),frequency,'facecolor',[0.6 0.6 0.6],'linewidth',0.01); +b=bar([1:max(bestModel)-1 max(bestModel)+1],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]; %meaning +b.CData(topModels(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 ylabel('fraction of cue cells'); -chance=1/90; +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 91]); +xlim([0 93]); ylim([0 0.3]); -xticks(topModels); +xticks(topModels(1:end-1)); yticks(0:0.1:0.3); %xtickangle(45); -%xticklabels(vnames(outliers)); +xticklabels(vnames(outliers(1:end-1))); +xtickangle(45); -category=8*ones(length(improvement),1); -category(cueK&~behavior)=7; %cue neurons +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 @@ -933,6 +969,7 @@ end end + %cueVector=cueVector-mean(cueVector,2); corrmat=corrcoef(cueVector); cueCorrMatrix(:,:,NN)=corrmat; @@ -1035,7 +1072,7 @@ end -titles={vnames{outliers},'non-meaning'}; +titles={vnames{outliers(1:end-1)},'untuned','odor'}; for ct=1:7 ax(2)=subplot(4,7,7+ct); @@ -1089,10 +1126,6 @@ end end -%% history coding - -[~,bestValue]=max(varExpValueSh(:,3:4),[],2); - %% remove doubles (from merging) and edge of FOV rois function iscell = processROIs(iscell,stat,F) ROIpixels={}; diff --git a/Imaging/dataInfo.txt b/Imaging/dataInfo.txt new file mode 100644 index 0000000..f849fe3 --- /dev/null +++ b/Imaging/dataInfo.txt @@ -0,0 +1 @@ +Imaging data available on figshare (see README) \ No newline at end of file diff --git a/Neuropixels/dataInfo.txt b/Neuropixels/dataInfo.txt new file mode 100644 index 0000000..b7c7c62 --- /dev/null +++ b/Neuropixels/dataInfo.txt @@ -0,0 +1 @@ +Electrophysiology data available on figshare (see README) \ No newline at end of file diff --git a/README.md b/README.md index 9b31f7a..a9b2b97 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ # ottenheimer-et-al-2022 Code and pointers to data associated with Ottenheimer, Hjort, Bowen, et al. 2022. -All code from preprint is uploaded. Data files will be uploaded separately. \ No newline at end of file +Code to analyze data is uploaded in this repository. Data files are available on figshare: https://doi.org/10.6084/m9.figshare.21365598 \ No newline at end of file