Browse Source

Upload files to 'AnalysisAndFigures'

David Ottenheimer 4 years ago
parent
commit
a5e57e2dcd

+ 345 - 0
AnalysisAndFigures/j_SingleUnitPavlovianAlcoholExposedAnalysis.m

@@ -0,0 +1,345 @@
+%%
+% Sept 26th, 2012: Added line 278. Now excludes Baseline calculations on excluded neurons
+% May 21st, 2012: Added the time of the max(PSTH) and the max for signigicant responses
+% May 15th, 2012: this version includes waveform analysis
+% May 9th, 2012: The baseline used to compute z-scores is matched with the prewindow matrix
+% May 9th, 2012: R.Erefnames added to avoid having to fetch it from the excel file
+% May 2nd, 2012: this version uses parametric tests (ttest/ttest2 instead of SignRank and Ranksum)
+% April 27th, 2012:
+%   -handles both GoNoGo and 2Go exps
+%   -uses the new cleaned-up ResDetectSignRank02 code
+% April 13th, 2012:
+%   -uses a different way of excluding neurons (it skips the excluded neurons instead of filtering then afterwards
+%   -R.Structure added that stores the brain region (DMS, CORE, SHELL)
+
+clear all; clc;
+global Dura Baseline Tm Tbase BSIZE Tbin
+tic
+path='RESULTpas.xls';
+load('RAWvppasALL.mat')
+RAW=RAWvppasALL;
+
+%Main settings
+SAVE_FLAG=1;
+BSIZE=0.01; %Do not change
+Dura=[-22 20]; Tm=Dura(1):BSIZE:Dura(2);
+%Baseline=[-22 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %now defined line 98
+Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
+PStat=0.05; %for comparing pre versus post windows, or event A versus event B
+MinNumTrials=1;
+RvppasALL=[];RvppasALL.Ninfo={};NN=0;Nneurons=0;
+BinSize=0.02;
+
+%Smoothing
+Smoothing=1; %0 for raw and 1 for smoothing
+SmoothTYPE='lowess';
+SmoothSPAN=5;
+if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end
+
+
+% List of events to analyze and analysis windows EXTRACTED from excel file
+[~,Erefnames]=xlsread(path,'Windows','a3:a12'); % cell that contains the event names
+prewin =  xlsread(path,'Windows','b3:c12');
+postwin = xlsread(path,'Windows','d3:e12');
+BLWin=    xlsread(path,'Windows','f3:g12');
+RespWin=  xlsread(path,'Windows','h3:i12');
+
+%Settings for Response detection w/ signrank
+WINin=[0.03,0.3]; %Onset and offset requirements in sec.
+%WINin=[-2 -4]; %negative values define onset requirement by the number of consecutive bins and not by duration
+CI=[.99,.99];
+
+%%
+%Finds the total number of neurons accross all sessions
+for i=1:length(RAW)
+    RvppasALL.Ninfo=cat(1,RvppasALL.Ninfo,RAW(i).Ninfo);
+    Nneurons=Nneurons+size(RAW(i).Nrast,1);
+end
+
+% load('CoordPPSall.mat');
+% RPPSall.Coord=CoordPPSall;
+RvppasALL.Coord=10*ones(Nneurons,4); %comment this line and uncomment the 2 lines above when you have real coordinates.
+RvppasALL.Structure=RvppasALL.Coord(:,4);
+RvppasALL.Ninfo=cat(2, RvppasALL.Ninfo, mat2cell([1:Nneurons]',ones(1,Nneurons),1));
+RvppasALL.Erefnames= Erefnames;
+
+RvppasALL.Rat(1:Nneurons,1)=NaN;
+RvppasALL.Session(1:Nneurons,1)=NaN;
+
+for i=1:length(RvppasALL.Ninfo)
+    RvppasALL.Rat(i,1)=str2num(RvppasALL.Ninfo{i,1}(4:5));
+    RvppasALL.Session(i,1)=str2num(RvppasALL.Ninfo{i,1}(10:11));
+end
+
+% preallocating
+RvppasALL.Param.Tm=Tm;
+RvppasALL.Param.Tbin=Tbin;
+RvppasALL.Param.Dura=Dura;
+RvppasALL.Param.Baseline=Baseline;
+RvppasALL.Param.PStat=PStat;
+RvppasALL.Param.MinNumTrials=MinNumTrials;
+RvppasALL.Param.path=path;
+RvppasALL.Param.prewin=prewin;
+RvppasALL.Param.postwin=postwin;
+RvppasALL.Param.BLwin=BLWin;
+RvppasALL.Param.RespWin=RespWin;
+RvppasALL.Param.ResponseReq=WINin;
+RvppasALL.Param.SmoothTYPE=SmoothTYPE;
+RvppasALL.Param.SmoothSPAN=SmoothSPAN;
+RvppasALL.CSPlusResponseStat(1:Nneurons,1)=NaN;
+RvppasALL.CSStat(1:Nneurons,1)=NaN;
+
+
+for k=1:length(Erefnames)
+    RvppasALL.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
+    RvppasALL.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
+    RvppasALL.Ev(k).Meanraw(1:Nneurons,1)=NaN;
+    RvppasALL.Ev(k).Meanz(1:Nneurons,1)=NaN;
+    RvppasALL.Ev(k).BW(1:Nneurons,1)=NaN;
+    RvppasALL.Ev(k).ttest(1:Nneurons,1)=NaN;
+    RvppasALL.Ev(k).RespDir(1:Nneurons,1)=NaN;
+    RvppasALL.Ev(k).RFLAG(1:Nneurons,1)=NaN;
+    RvppasALL.Ev(k).CFLAG(1:Nneurons,1)=NaN;
+    RvppasALL.Ev(k).Onsets(1:Nneurons,:)=NaN(Nneurons,2);
+    RvppasALL.Ev(k).Offsets(1:Nneurons,:)=NaN(Nneurons,2);
+    RvppasALL.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
+    RvppasALL.Ev(k).MaxVal(1:Nneurons,1)=NaN;
+    RvppasALL.Ev(k).MaxTime(1:Nneurons,1)=NaN;
+end
+
+%% runs the main routine
+
+for i=1:length(RAW) %loops through sessions
+    for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
+        NN=NN+1; %neuron counter
+        if RvppasALL.Structure(NN)~=0 %to avoid analyzing excluded neurons
+            for k=1:length(Erefnames) %loops thorough the events
+                EvInd=strcmp(Erefnames(k),RAW(i).Einfo(:,2)); %find the event id number from RAW
+                
+                if sum(EvInd)==0
+                    fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
+                end
+                
+                RvppasALL.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
+                Tbase=prewin(k,1):BSIZE:prewin(k,2);
+                if  ~isempty(EvInd) && RvppasALL.Ev(k).NumberTrials(NN,1)>=MinNumTrials %avoid analyzing sessions where that do not have enough trials
+                    [PSR0,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{1});% makes collapsed rasters for baseline.
+                    [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
+                    [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
+                    if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
+                        %optimal bin size
+                        %[PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,1});%-----DP used here
+                        %[PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BW1});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
+                        %Fixed bin size
+                        [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,0,BinSize});%-----DP used here
+                        [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
+                        PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
+                        PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';
+                        
+                        %------------- Fills the R.Ev(k) fields --------------
+                        RvppasALL.Ev(k).BW(NN,1)=BW1;
+                        RvppasALL.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
+                        RvppasALL.Ev(k).Meanraw(NN,1)=nanmean(RvppasALL.Ev(k).PSTHraw(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
+                        if sum(PTH0,2)~=0
+                            RvppasALL.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
+                            RvppasALL.Ev(k).Meanz(NN,1)=nanmean(RvppasALL.Ev(k).PSTHz(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
+                        else
+                            RvppasALL.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
+                            RvppasALL.Ev(k).Meanz(NN,1)=NaN;
+                        end
+                        
+                        %------------------ firing (in Hz) per trial in pre/post windows ------------------
+                        %used to make the between events comparisons and Response detection in a single window----
+                        ev(k).pre=NaN(size(RAW(i).Erast{EvInd},1),1);
+                        ev(k).post=NaN(size(RAW(i).Erast{EvInd},1),1);
+                        for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
+                            ev(k).pre(m)=sum(PSR2{m}<prewin(k,2) & PSR2{m}>prewin(k,1))/(prewin(k,2)-prewin(k,1));
+                            ev(k).post(m)=sum(PSR2{m}<postwin(k,2) & PSR2{m}>postwin(k,1))/(postwin(k,2)-postwin(k,1));
+                        end
+                        
+                        %---------------------------Response detection w/ SignRank on pre/post windows---------------------------
+                        [~,RvppasALL.Ev(k).ttest(NN,1)]=ttest(ev(k).pre, ev(k).post); %Signrank used here because it is a dependant sample test
+                        if RvppasALL.Ev(k).ttest(NN,1)<PStat
+                            RvppasALL.Ev(k).RespDir(NN,1)=sign(mean(ev(k).post)-mean(ev(k).pre));
+                            if RvppasALL.Ev(k).RespDir(NN,1)>0
+                                Search=find(Tm>=postwin(k,1) & Tm<=postwin(k,2));
+                                [RvppasALL.Ev(k).MaxVal(NN,1),MaxInd]=max(RvppasALL.Ev(k).PSTHraw(NN,Search));
+                                RvppasALL.Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
+                            else
+                                Search=find(Tm>=postwin(k,1) & Tm<=postwin(k,2));
+                                [RvppasALL.Ev(k).MaxVal(NN,1),MaxInd]=min(RvppasALL.Ev(k).PSTHraw(NN,Search));
+                                RvppasALL.Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
+                            end
+                        else RvppasALL.Ev(k).RespDir(NN,1)=0;
+                        end
+                        
+                        %---------------------------Response detection w/ RESDETECT08 on running windows---------------------------
+                        RD=ResDetect08(PTH1,PTH0,RespWin(k,:),0,WINin, CI);
+                        
+                        RvppasALL.Ev(k).RFLAG(NN,1)=RD.RFLAG;
+                        RvppasALL.Ev(k).CFLAG(NN,1)=RD.CFLAG;
+                        RvppasALL.Ev(k).Onsets(NN,:)=RD.Onsets';
+                        RvppasALL.Ev(k).Offsets(NN,:)=RD.Offsets';
+                        RvppasALL.Param.Paramnames=RD.Paramnames;
+                        RvppasALL.Param.RDParam=RD.Params;
+                        
+                    end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
+                end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
+            end %Events
+            
+            
+            %----------------------- CUES -----------------------
+            CSPlus=strcmp('CSPlus', Erefnames);
+            CSMinus=strcmp('CSMinus', Erefnames);
+            if sum(CSPlus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+                if sum(CSMinus)~=0
+                    RvppasALL.CScriteria=linspace(0, max([ev(CSPlus).post(:);ev(CSMinus).post(:)]),200);%Determine the criteria range
+                    RvppasALL.CSprecriteria=linspace(0, max([ev(CSPlus).pre(:);ev(CSMinus).pre(:)]),200);%Determine the criteria range
+                    for s=1:length(RvppasALL.CScriteria);
+                        RvppasALL.CSfalsepos(NN,s)=sum(ev(CSMinus).post>=RvppasALL.CScriteria(s))/length(ev(CSMinus).post); %Determine the probability that the baseline is greater than each criteria (x)
+                        RvppasALL.CStruepos(NN,s)=sum(ev(CSPlus).post>=RvppasALL.CScriteria(s))/length(ev(CSPlus).post);%Determine the probability that the post-window firing firing is greater than critera(y)
+                        RvppasALL.CSprefalsepos(NN,s)=sum(ev(CSMinus).pre>=RvppasALL.CScriteria(s))/length(ev(CSMinus).pre); %Determine the probability that the baseline is greater than each criteria (x)
+                        RvppasALL.CSpretruepos(NN,s)=sum(ev(CSPlus).pre>=RvppasALL.CScriteria(s))/length(ev(CSPlus).pre);%Determine the probability that the pre-window firing firing is greater than critera(y)
+                    end
+                    RvppasALL.CSauROC(NN,1)=trapz(RvppasALL.CSfalsepos(NN,:),RvppasALL.CStruepos(NN,:));%Calculate area under the curve for this neuron for this event
+                    RvppasALL.CSpreauROC(NN,1)=trapz(RvppasALL.CSprefalsepos(NN,:),RvppasALL.CSpretruepos(NN,:));%Calculate area under the curve for this neuron for this event
+                    [~,RvppasALL.CSStat(NN,1)]=ttest2(ev(CSPlus).post,ev(CSMinus).post); %Ranksum test used becasue it is an independant sample test
+                else RvppasALL.CSStat(NN,1)=NaN;
+                     RvppasALL.CSauROC(NN,1)=NaN;
+                     RvppasALL.CSpreauROC(NN,1)=NaN;
+                end
+            end
+ 
+                                
+            
+            PECSPlus=strcmp('PECSPlus', Erefnames);
+            PECSMinus=strcmp('PECSMinus', Erefnames);
+            if sum(PECSPlus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+                if sum(PECSMinus)~=0
+                    [RvppasALL.PECueStat(NN,1),~]=ranksum(ev(PECSPlus).post,ev(PECSMinus).post); %Ranksum test used becasue it is an independant sample test
+                else RvppasALL.PECueStat(NN,1)=NaN;
+                end
+            end
+            
+            PECSPlus=strcmp('PECSPlus', Erefnames);
+            ITIPE=strcmp('ITIPE', Erefnames);
+            if sum(PECSPlus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+                if RvppasALL.Ev(PECSPlus).RespDir(NN,1)~=0 || RvppasALL.Ev(ITIPE).RespDir(NN,1)~=0
+                    [RvppasALL.PERewStat(NN,1),~]=ranksum(ev(PECSPlus).post,ev(ITIPE).post); %Ranksum test used becasue it is an independant sample test
+                else RvppasALL.PERewStat(NN,1)=NaN;
+                end
+            end
+            
+            %             PECSMinus=strcmp('PECSMinus', Erefnames);
+            %             ITIPE=strcmp('ITIPE', Erefnames);
+            %             if sum(PECSMinus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+            %                 if R.Ev(PECSMinus).RespDir(NN,1)~=0 || R.Ev(ITIPE).RespDir(NN,1)~=0
+            %                     [R.PEUnrewCueStat(NN,1),~]=ranksum(ev(PECSMinus).post,ev(ITIPE).post); %Ranksum test used becasue it is an independant sample test
+            %                 else R.PEUnrewCueStat(NN,1)=NaN;
+            %                 end
+            %             end
+            
+            CSPluswPE=strcmp('CSPluswPE', Erefnames);
+            CSPlusnoPE=strcmp('CSPlusnoPE', Erefnames);
+            if sum(CSPlusnoPE)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+                if sum(CSPluswPE)~=0
+                    RvppasALL.CSrespcriteria=linspace(0, max([ev(CSPlusnoPE).post(:);ev(CSPluswPE).post(:)]),200);%Determine the criteria range
+                    RvppasALL.CSprerespcriteria=linspace(0, max([ev(CSPlusnoPE).pre(:);ev(CSPluswPE).pre(:)]),200);%Determine the criteria range
+                    for s=1:length(RvppasALL.CSrespcriteria);
+                        RvppasALL.CSrespfalsepos(NN,s)=sum(ev(CSPlusnoPE).post>=RvppasALL.CSrespcriteria(s))/length(ev(CSPlusnoPE).post); %Determine the probability that the baseline is greater than each criteria (x)
+                        RvppasALL.CSresptruepos(NN,s)=sum(ev(CSPluswPE).post>=RvppasALL.CSrespcriteria(s))/length(ev(CSPluswPE).post);%Determine the probability that the post-window firing firing is greater than critera(y)
+                        RvppasALL.CSprerespfalsepos(NN,s)=sum(ev(CSPlusnoPE).pre>=RvppasALL.CSrespcriteria(s))/length(ev(CSPlusnoPE).pre); %Determine the probability that the baseline is greater than each criteria (x)
+                        RvppasALL.CSpreresptruepos(NN,s)=sum(ev(CSPluswPE).pre>=RvppasALL.CSrespcriteria(s))/length(ev(CSPluswPE).pre);%Determine the probability that the post-window firing firing is greater than critera(y)
+                    end
+                    RvppasALL.CSrespauROC(NN,1)=trapz(RvppasALL.CSrespfalsepos(NN,:),RvppasALL.CSresptruepos(NN,:));%Calculate area under the curve for this neuron for this event
+                    RvppasALL.CSprerespauROC(NN,1)=trapz(RvppasALL.CSprerespfalsepos(NN,:),RvppasALL.CSpreresptruepos(NN,:));
+                    [~,RvppasALL.CSPlusResponseStat(NN,1)]=ttest2(ev(CSPluswPE).post,ev(CSPlusnoPE).post); %Ranksum test used becasue it is an independant sample test
+                else RvppasALL.CSPlusResponseStat(NN,1)=NaN;
+                    RvppasALL.CSrespauROC(NN,1)=NaN;
+                end
+            else
+                RvppasALL.CSrespauROC(NN,1)=NaN;
+                RvppasALL.CSPlusResponseStat(NN,1)=NaN;
+                RvppasALL.CSprerespauROC(NN,1)=NaN;
+            end
+            
+            
+            
+            %
+            CSMinuswPE=strcmp('CSMinuswPE', Erefnames);
+            CSMinusnoPE=strcmp('CSMinusnoPE', Erefnames);
+            if sum(CSMinuswPE)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+                if RvppasALL.Ev(CSMinuswPE).RespDir(NN,1)~=0 || RvppasALL.Ev(CSMinusnoPE).RespDir(NN,1)~=0
+                    [RvppasALL.CSMinusResponseStat(NN,1),~]=ranksum(ev(CSMinuswPE).post,ev(CSMinusnoPE).post); %Ranksum test used becasue it is an independant sample test
+                else RvppasALL.CSMinusResponseStat(NN,1)=NaN;
+                end
+            end
+            %
+    RvppasALL.CSPlusRatio(NN,1)=length(ev(CSPluswPE).post)/length(ev(CSPlus).post);
+    RvppasALL.CSMinusRatio(NN,1)=length(ev(CSMinuswPE).post)/length(ev(CSMinus).post);
+            fprintf('Neuron ID # %d\n',NN);
+        elseif RvppasALL.Structure(NN)~=0
+        end %exclusion: IF R.Structure(NN)~=0 to avoid analyzing excluded neurons
+    end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
+end %sessions: FOR i=1:length(RAW)
+
+RvppasALL.CSfalseposMEAN(1,:)=nanmean(RvppasALL.CSfalsepos(RvppasALL.Structure==10,:));
+RvppasALL.CSfalseposSEM(1,:)=nanste(RvppasALL.CSfalsepos(RvppasALL.Structure==10,:),1);
+
+RvppasALL.CStrueposMEAN(1,:)=nanmean(RvppasALL.CStruepos(RvppasALL.Structure==10,:));
+RvppasALL.CStrueposSEM(1,:)=nanste(RvppasALL.CStruepos(RvppasALL.Structure==10,:),1);
+
+RvppasALL.CSrespfalseposMEAN(1,:)=nanmean(RvppasALL.CSrespfalsepos(RvppasALL.Structure==10,:));
+RvppasALL.CSrespfalseposSEM(1,:)=nanste(RvppasALL.CSrespfalsepos(RvppasALL.Structure==10,:),1);
+
+RvppasALL.CSresptrueposMEAN(1,:)=nanmean(RvppasALL.CSresptruepos(RvppasALL.Structure==10,:));
+RvppasALL.CSresptrueposSEM(1,:)=nanste(RvppasALL.CSresptruepos(RvppasALL.Structure==10,:),1);
+
+% RPPSall.CSprefalseposMEAN(1,:)=nanmean(RPPSall.CSprefalsepos(RPPSall.Structure==10,:));
+% RPPSall.CSprefalseposSEM(1,:)=nanste(RPPSall.CSprefalsepos(RPPSall.Structure==10,:),1);
+% 
+% RPPSall.CSpretrueposMEAN(1,:)=nanmean(RPPSall.CSpretruepos(RPPSall.Structure==10,:));
+% RPPSall.CSpretrueposSEM(1,:)=nanste(RPPSall.CSpretruepos(RPPSall.Structure==10,:),1);
+
+
+if SAVE_FLAG
+    save('RvppasALL.mat','RvppasALL')
+end
+
+toc
+
+
+%% Create a Summury excel spreadsheet
+
+% RESULT=[];
+% %RESULT=cat(2,RESULT, R.CueCorStat, R.CueCorRExtStat,R.CueIncorStat,R.CueIncorRExtStat,R.BaselineTrend, R.LExtTrend, R.BaselineLExtOne, R.BaselineLExtNo, R.DeltaLP, R.DeltaNo);
+% 
+% for k=1:length(Erefnames)
+%     RESULT=cat(2,RESULT,RPPSall.Ev(k).RespDir, RPPSall.Ev(k).RFLAG, RPPSall.Ev(k).CFLAG,RPPSall.Ev(k).Onsets, RPPSall.Ev(k).Offsets);
+% end
+% 
+% xlswrite(path,RPPSall.Ninfo,'RESULTsheet', 'A3'); % Writes Session - Neuron - ID#
+% xlswrite(path,RPPSall.Structure,'RESULTsheet', 'D3');  % Writes the brain region
+% xlswrite(path,RESULT,'RESULTsheet', 'F3'); % Writes RESULT matrix
+% toc
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

+ 345 - 0
AnalysisAndFigures/k_SingleUnitInstrumentalAlcoholExposedAnalysis.m

@@ -0,0 +1,345 @@
+%%
+% Sept 26th, 2012: Added line 278. Now excludes Baseline calculations on excluded neurons
+% May 21st, 2012: Added the time of the max(PSTH) and the max for signigicant responses
+% May 15th, 2012: this version includes waveform analysis
+% May 9th, 2012: The baseline used to compute z-scores is matched with the prewindow matrix
+% May 9th, 2012: R.Erefnames added to avoid having to fetch it from the excel file
+% May 2nd, 2012: this version uses parametric tests (ttest/ttest2 instead of SignRank and Ranksum)
+% April 27th, 2012:
+%   -handles both GoNoGo and 2Go exps
+%   -uses the new cleaned-up ResDetectSignRank02 code
+% April 13th, 2012:
+%   -uses a different way of excluding neurons (it skips the excluded neurons instead of filtering then afterwards
+%   -R.Structure added that stores the brain region (DMS, CORE, SHELL)
+
+clear all; clc;
+global Dura Baseline Tm Tbase BSIZE Tbin
+tic
+path='RESULTppas.xls';
+load('RAWvpppasALL.mat')
+RAW=RAWvpppasALL;
+
+%Main settings
+SAVE_FLAG=1;
+BSIZE=0.01; %Do not change
+Dura=[-22 20]; Tm=Dura(1):BSIZE:Dura(2);
+%Baseline=[-22 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %now defined line 98
+Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
+PStat=0.05; %for comparing pre versus post windows, or event A versus event B
+MinNumTrials=1;
+RvpppasALL=[];RvpppasALL.Ninfo={};NN=0;Nneurons=0;
+BinSize=0.02;
+
+%Smoothing
+Smoothing=1; %0 for raw and 1 for smoothing
+SmoothTYPE='lowess';
+SmoothSPAN=5;
+if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end
+
+
+% List of events to analyze and analysis windows EXTRACTED from excel file
+[~,Erefnames]=xlsread(path,'Windows','a3:a12'); % cell that contains the event names
+prewin =  xlsread(path,'Windows','b3:c12');
+postwin = xlsread(path,'Windows','d3:e12');
+BLWin=    xlsread(path,'Windows','f3:g12');
+RespWin=  xlsread(path,'Windows','h3:i12');
+
+%Settings for Response detection w/ signrank
+WINin=[0.03,0.3]; %Onset and offset requirements in sec.
+%WINin=[-2 -4]; %negative values define onset requirement by the number of consecutive bins and not by duration
+CI=[.99,.99];
+
+%%
+%Finds the total number of neurons accross all sessions
+for i=1:length(RAW)
+    RvpppasALL.Ninfo=cat(1,RvpppasALL.Ninfo,RAW(i).Ninfo);
+    Nneurons=Nneurons+size(RAW(i).Nrast,1);
+end
+
+% load('CoordPPSall.mat');
+% RPPSall.Coord=CoordPPSall;
+RvpppasALL.Coord=10*ones(Nneurons,4); %comment this line and uncomment the 2 lines above when you have real coordinates.
+RvpppasALL.Structure=RvpppasALL.Coord(:,4);
+RvpppasALL.Ninfo=cat(2, RvpppasALL.Ninfo, mat2cell([1:Nneurons]',ones(1,Nneurons),1));
+RvpppasALL.Erefnames= Erefnames;
+
+RvpppasALL.Rat(1:Nneurons,1)=NaN;
+RvpppasALL.Session(1:Nneurons,1)=NaN;
+
+for i=1:length(RvpppasALL.Ninfo)
+    RvpppasALL.Rat(i,1)=str2num(RvpppasALL.Ninfo{i,1}(4:5));
+    RvpppasALL.Session(i,1)=str2num(RvpppasALL.Ninfo{i,1}(10:11));
+end
+
+% preallocating
+RvpppasALL.Param.Tm=Tm;
+RvpppasALL.Param.Tbin=Tbin;
+RvpppasALL.Param.Dura=Dura;
+RvpppasALL.Param.Baseline=Baseline;
+RvpppasALL.Param.PStat=PStat;
+RvpppasALL.Param.MinNumTrials=MinNumTrials;
+RvpppasALL.Param.path=path;
+RvpppasALL.Param.prewin=prewin;
+RvpppasALL.Param.postwin=postwin;
+RvpppasALL.Param.BLwin=BLWin;
+RvpppasALL.Param.RespWin=RespWin;
+RvpppasALL.Param.ResponseReq=WINin;
+RvpppasALL.Param.SmoothTYPE=SmoothTYPE;
+RvpppasALL.Param.SmoothSPAN=SmoothSPAN;
+RvpppasALL.CSPlusResponseStat(1:Nneurons,1)=NaN;
+RvpppasALL.CSStat(1:Nneurons,1)=NaN;
+
+
+for k=1:length(Erefnames)
+    RvpppasALL.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
+    RvpppasALL.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
+    RvpppasALL.Ev(k).Meanraw(1:Nneurons,1)=NaN;
+    RvpppasALL.Ev(k).Meanz(1:Nneurons,1)=NaN;
+    RvpppasALL.Ev(k).BW(1:Nneurons,1)=NaN;
+    RvpppasALL.Ev(k).ttest(1:Nneurons,1)=NaN;
+    RvpppasALL.Ev(k).RespDir(1:Nneurons,1)=NaN;
+    RvpppasALL.Ev(k).RFLAG(1:Nneurons,1)=NaN;
+    RvpppasALL.Ev(k).CFLAG(1:Nneurons,1)=NaN;
+    RvpppasALL.Ev(k).Onsets(1:Nneurons,:)=NaN(Nneurons,2);
+    RvpppasALL.Ev(k).Offsets(1:Nneurons,:)=NaN(Nneurons,2);
+    RvpppasALL.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
+    RvpppasALL.Ev(k).MaxVal(1:Nneurons,1)=NaN;
+    RvpppasALL.Ev(k).MaxTime(1:Nneurons,1)=NaN;
+end
+
+%% runs the main routine
+
+for i=1:length(RAW) %loops through sessions
+    for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
+        NN=NN+1; %neuron counter
+        if RvpppasALL.Structure(NN)~=0 %to avoid analyzing excluded neurons
+            for k=1:length(Erefnames) %loops thorough the events
+                EvInd=strcmp(Erefnames(k),RAW(i).Einfo(:,2)); %find the event id number from RAW
+                
+                if sum(EvInd)==0
+                    fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
+                end
+                
+                RvpppasALL.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
+                Tbase=prewin(k,1):BSIZE:prewin(k,2);
+                if  ~isempty(EvInd) && RvpppasALL.Ev(k).NumberTrials(NN,1)>=MinNumTrials %avoid analyzing sessions where that do not have enough trials
+                    [PSR0,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{1});% makes collapsed rasters for baseline.
+                    [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
+                    [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
+                    if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
+                        %optimal bin size
+                        %[PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,1});%-----DP used here
+                        %[PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BW1});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
+                        %Fixed bin size
+                        [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,0,BinSize});%-----DP used here
+                        [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
+                        PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
+                        PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';
+                        
+                        %------------- Fills the R.Ev(k) fields --------------
+                        RvpppasALL.Ev(k).BW(NN,1)=BW1;
+                        RvpppasALL.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
+                        RvpppasALL.Ev(k).Meanraw(NN,1)=nanmean(RvpppasALL.Ev(k).PSTHraw(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
+                        if sum(PTH0,2)~=0
+                            RvpppasALL.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
+                            RvpppasALL.Ev(k).Meanz(NN,1)=nanmean(RvpppasALL.Ev(k).PSTHz(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
+                        else
+                            RvpppasALL.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
+                            RvpppasALL.Ev(k).Meanz(NN,1)=NaN;
+                        end
+                        
+                        %------------------ firing (in Hz) per trial in pre/post windows ------------------
+                        %used to make the between events comparisons and Response detection in a single window----
+                        ev(k).pre=NaN(size(RAW(i).Erast{EvInd},1),1);
+                        ev(k).post=NaN(size(RAW(i).Erast{EvInd},1),1);
+                        for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
+                            ev(k).pre(m)=sum(PSR2{m}<prewin(k,2) & PSR2{m}>prewin(k,1))/(prewin(k,2)-prewin(k,1));
+                            ev(k).post(m)=sum(PSR2{m}<postwin(k,2) & PSR2{m}>postwin(k,1))/(postwin(k,2)-postwin(k,1));
+                        end
+                        
+                        %---------------------------Response detection w/ SignRank on pre/post windows---------------------------
+                        [~,RvpppasALL.Ev(k).ttest(NN,1)]=ttest(ev(k).pre, ev(k).post); %Signrank used here because it is a dependant sample test
+                        if RvpppasALL.Ev(k).ttest(NN,1)<PStat
+                            RvpppasALL.Ev(k).RespDir(NN,1)=sign(mean(ev(k).post)-mean(ev(k).pre));
+                            if RvpppasALL.Ev(k).RespDir(NN,1)>0
+                                Search=find(Tm>=postwin(k,1) & Tm<=postwin(k,2));
+                                [RvpppasALL.Ev(k).MaxVal(NN,1),MaxInd]=max(RvpppasALL.Ev(k).PSTHraw(NN,Search));
+                                RvpppasALL.Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
+                            else
+                                Search=find(Tm>=postwin(k,1) & Tm<=postwin(k,2));
+                                [RvpppasALL.Ev(k).MaxVal(NN,1),MaxInd]=min(RvpppasALL.Ev(k).PSTHraw(NN,Search));
+                                RvpppasALL.Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
+                            end
+                        else RvpppasALL.Ev(k).RespDir(NN,1)=0;
+                        end
+                        
+                        %---------------------------Response detection w/ RESDETECT08 on running windows---------------------------
+                        RD=ResDetect08(PTH1,PTH0,RespWin(k,:),0,WINin, CI);
+                        
+                        RvpppasALL.Ev(k).RFLAG(NN,1)=RD.RFLAG;
+                        RvpppasALL.Ev(k).CFLAG(NN,1)=RD.CFLAG;
+                        RvpppasALL.Ev(k).Onsets(NN,:)=RD.Onsets';
+                        RvpppasALL.Ev(k).Offsets(NN,:)=RD.Offsets';
+                        RvpppasALL.Param.Paramnames=RD.Paramnames;
+                        RvpppasALL.Param.RDParam=RD.Params;
+                        
+                    end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
+                end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
+            end %Events
+            
+            
+            %----------------------- CUES -----------------------
+            CSPlus=strcmp('CSPlus', Erefnames);
+            CSMinus=strcmp('CSMinus', Erefnames);
+            if sum(CSPlus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+                if sum(CSMinus)~=0
+                    RvpppasALL.CScriteria=linspace(0, max([ev(CSPlus).post(:);ev(CSMinus).post(:)]),200);%Determine the criteria range
+                    RvpppasALL.CSprecriteria=linspace(0, max([ev(CSPlus).pre(:);ev(CSMinus).pre(:)]),200);%Determine the criteria range
+                    for s=1:length(RvpppasALL.CScriteria);
+                        RvpppasALL.CSfalsepos(NN,s)=sum(ev(CSMinus).post>=RvpppasALL.CScriteria(s))/length(ev(CSMinus).post); %Determine the probability that the baseline is greater than each criteria (x)
+                        RvpppasALL.CStruepos(NN,s)=sum(ev(CSPlus).post>=RvpppasALL.CScriteria(s))/length(ev(CSPlus).post);%Determine the probability that the post-window firing firing is greater than critera(y)
+                        RvpppasALL.CSprefalsepos(NN,s)=sum(ev(CSMinus).pre>=RvpppasALL.CScriteria(s))/length(ev(CSMinus).pre); %Determine the probability that the baseline is greater than each criteria (x)
+                        RvpppasALL.CSpretruepos(NN,s)=sum(ev(CSPlus).pre>=RvpppasALL.CScriteria(s))/length(ev(CSPlus).pre);%Determine the probability that the pre-window firing firing is greater than critera(y)
+                    end
+                    RvpppasALL.CSauROC(NN,1)=trapz(RvpppasALL.CSfalsepos(NN,:),RvpppasALL.CStruepos(NN,:));%Calculate area under the curve for this neuron for this event
+                    RvpppasALL.CSpreauROC(NN,1)=trapz(RvpppasALL.CSprefalsepos(NN,:),RvpppasALL.CSpretruepos(NN,:));%Calculate area under the curve for this neuron for this event
+                    [~,RvpppasALL.CSStat(NN,1)]=ttest2(ev(CSPlus).post,ev(CSMinus).post); %Ranksum test used becasue it is an independant sample test
+                else RvpppasALL.CSStat(NN,1)=NaN;
+                     RvpppasALL.CSauROC(NN,1)=NaN;
+                     RvpppasALL.CSpreauROC(NN,1)=NaN;
+                end
+            end
+ 
+                                
+            
+            PECSPlus=strcmp('PECSPlus', Erefnames);
+            PECSMinus=strcmp('PECSMinus', Erefnames);
+            if sum(PECSPlus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+                if sum(PECSMinus)~=0
+                    [RvpppasALL.PECueStat(NN,1),~]=ranksum(ev(PECSPlus).post,ev(PECSMinus).post); %Ranksum test used becasue it is an independant sample test
+                else RvpppasALL.PECueStat(NN,1)=NaN;
+                end
+            end
+            
+            PECSPlus=strcmp('PECSPlus', Erefnames);
+            ITIPE=strcmp('ITIPE', Erefnames);
+            if sum(PECSPlus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+                if RvpppasALL.Ev(PECSPlus).RespDir(NN,1)~=0 || RvpppasALL.Ev(ITIPE).RespDir(NN,1)~=0
+                    [RvpppasALL.PERewStat(NN,1),~]=ranksum(ev(PECSPlus).post,ev(ITIPE).post); %Ranksum test used becasue it is an independant sample test
+                else RvpppasALL.PERewStat(NN,1)=NaN;
+                end
+            end
+            
+            %             PECSMinus=strcmp('PECSMinus', Erefnames);
+            %             ITIPE=strcmp('ITIPE', Erefnames);
+            %             if sum(PECSMinus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+            %                 if R.Ev(PECSMinus).RespDir(NN,1)~=0 || R.Ev(ITIPE).RespDir(NN,1)~=0
+            %                     [R.PEUnrewCueStat(NN,1),~]=ranksum(ev(PECSMinus).post,ev(ITIPE).post); %Ranksum test used becasue it is an independant sample test
+            %                 else R.PEUnrewCueStat(NN,1)=NaN;
+            %                 end
+            %             end
+            
+            CSPluswPE=strcmp('CSPluswPE', Erefnames);
+            CSPlusnoPE=strcmp('CSPlusnoPE', Erefnames);
+            if sum(CSPlusnoPE)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+                if sum(CSPluswPE)~=0
+                    RvpppasALL.CSrespcriteria=linspace(0, max([ev(CSPlusnoPE).post(:);ev(CSPluswPE).post(:)]),200);%Determine the criteria range
+                    RvpppasALL.CSprerespcriteria=linspace(0, max([ev(CSPlusnoPE).pre(:);ev(CSPluswPE).pre(:)]),200);%Determine the criteria range
+                    for s=1:length(RvpppasALL.CSrespcriteria);
+                        RvpppasALL.CSrespfalsepos(NN,s)=sum(ev(CSPlusnoPE).post>=RvpppasALL.CSrespcriteria(s))/length(ev(CSPlusnoPE).post); %Determine the probability that the baseline is greater than each criteria (x)
+                        RvpppasALL.CSresptruepos(NN,s)=sum(ev(CSPluswPE).post>=RvpppasALL.CSrespcriteria(s))/length(ev(CSPluswPE).post);%Determine the probability that the post-window firing firing is greater than critera(y)
+                        RvpppasALL.CSprerespfalsepos(NN,s)=sum(ev(CSPlusnoPE).pre>=RvpppasALL.CSrespcriteria(s))/length(ev(CSPlusnoPE).pre); %Determine the probability that the baseline is greater than each criteria (x)
+                        RvpppasALL.CSpreresptruepos(NN,s)=sum(ev(CSPluswPE).pre>=RvpppasALL.CSrespcriteria(s))/length(ev(CSPluswPE).pre);%Determine the probability that the post-window firing firing is greater than critera(y)
+                    end
+                    RvpppasALL.CSrespauROC(NN,1)=trapz(RvpppasALL.CSrespfalsepos(NN,:),RvpppasALL.CSresptruepos(NN,:));%Calculate area under the curve for this neuron for this event
+                    RvpppasALL.CSprerespauROC(NN,1)=trapz(RvpppasALL.CSprerespfalsepos(NN,:),RvpppasALL.CSpreresptruepos(NN,:));
+                    [~,RvpppasALL.CSPlusResponseStat(NN,1)]=ttest2(ev(CSPluswPE).post,ev(CSPlusnoPE).post); %Ranksum test used becasue it is an independant sample test
+                else RvpppasALL.CSPlusResponseStat(NN,1)=NaN;
+                    RvpppasALL.CSrespauROC(NN,1)=NaN;
+                end
+            else
+                RvpppasALL.CSrespauROC(NN,1)=NaN;
+                RvpppasALL.CSPlusResponseStat(NN,1)=NaN;
+                RvpppasALL.CSprerespauROC(NN,1)=NaN;
+            end
+            
+            
+            
+            %
+            CSMinuswPE=strcmp('CSMinuswPE', Erefnames);
+            CSMinusnoPE=strcmp('CSMinusnoPE', Erefnames);
+            if sum(CSMinuswPE)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
+                if RvpppasALL.Ev(CSMinuswPE).RespDir(NN,1)~=0 || RvpppasALL.Ev(CSMinusnoPE).RespDir(NN,1)~=0
+                    [RvpppasALL.CSMinusResponseStat(NN,1),~]=ranksum(ev(CSMinuswPE).post,ev(CSMinusnoPE).post); %Ranksum test used becasue it is an independant sample test
+                else RvpppasALL.CSMinusResponseStat(NN,1)=NaN;
+                end
+            end
+            %
+    RvpppasALL.CSPlusRatio(NN,1)=length(ev(CSPluswPE).post)/length(ev(CSPlus).post);
+    RvpppasALL.CSMinusRatio(NN,1)=length(ev(CSMinuswPE).post)/length(ev(CSMinus).post);
+            fprintf('Neuron ID # %d\n',NN);
+        elseif RvpppasALL.Structure(NN)~=0
+        end %exclusion: IF R.Structure(NN)~=0 to avoid analyzing excluded neurons
+    end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
+end %sessions: FOR i=1:length(RAW)
+
+RvpppasALL.CSfalseposMEAN(1,:)=nanmean(RvpppasALL.CSfalsepos(RvpppasALL.Structure==10,:));
+RvpppasALL.CSfalseposSEM(1,:)=nanste(RvpppasALL.CSfalsepos(RvpppasALL.Structure==10,:),1);
+
+RvpppasALL.CStrueposMEAN(1,:)=nanmean(RvpppasALL.CStruepos(RvpppasALL.Structure==10,:));
+RvpppasALL.CStrueposSEM(1,:)=nanste(RvpppasALL.CStruepos(RvpppasALL.Structure==10,:),1);
+
+RvpppasALL.CSrespfalseposMEAN(1,:)=nanmean(RvpppasALL.CSrespfalsepos(RvpppasALL.Structure==10,:));
+RvpppasALL.CSrespfalseposSEM(1,:)=nanste(RvpppasALL.CSrespfalsepos(RvpppasALL.Structure==10,:),1);
+
+RvpppasALL.CSresptrueposMEAN(1,:)=nanmean(RvpppasALL.CSresptruepos(RvpppasALL.Structure==10,:));
+RvpppasALL.CSresptrueposSEM(1,:)=nanste(RvpppasALL.CSresptruepos(RvpppasALL.Structure==10,:),1);
+
+% RPPSall.CSprefalseposMEAN(1,:)=nanmean(RPPSall.CSprefalsepos(RPPSall.Structure==10,:));
+% RPPSall.CSprefalseposSEM(1,:)=nanste(RPPSall.CSprefalsepos(RPPSall.Structure==10,:),1);
+% 
+% RPPSall.CSpretrueposMEAN(1,:)=nanmean(RPPSall.CSpretruepos(RPPSall.Structure==10,:));
+% RPPSall.CSpretrueposSEM(1,:)=nanste(RPPSall.CSpretruepos(RPPSall.Structure==10,:),1);
+
+
+if SAVE_FLAG
+    save('RvpppasALL.mat','RPPASall')
+end
+
+toc
+
+
+%% Create a Summury excel spreadsheet
+
+% RESULT=[];
+% %RESULT=cat(2,RESULT, R.CueCorStat, R.CueCorRExtStat,R.CueIncorStat,R.CueIncorRExtStat,R.BaselineTrend, R.LExtTrend, R.BaselineLExtOne, R.BaselineLExtNo, R.DeltaLP, R.DeltaNo);
+% 
+% for k=1:length(Erefnames)
+%     RESULT=cat(2,RESULT,RPPSall.Ev(k).RespDir, RPPSall.Ev(k).RFLAG, RPPSall.Ev(k).CFLAG,RPPSall.Ev(k).Onsets, RPPSall.Ev(k).Offsets);
+% end
+% 
+% xlswrite(path,RPPSall.Ninfo,'RESULTsheet', 'A3'); % Writes Session - Neuron - ID#
+% xlswrite(path,RPPSall.Structure,'RESULTsheet', 'D3');  % Writes the brain region
+% xlswrite(path,RESULT,'RESULTsheet', 'F3'); % Writes RESULT matrix
+% toc
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

+ 668 - 0
AnalysisAndFigures/l_Figure5_Figure6_BehaviorAndSingleUnitActivityAlcoholExposed.m

@@ -0,0 +1,668 @@
+%% Code for final figures for pavlovian versus instrumental comparison
+clear all;
+if ~exist('RvppasALL'),load('RvppasALL.mat'); end
+if ~exist('RvpppasALL'),load('RvpppasALL.mat'); end
+DarkPurple=[123/255 50/255 148/255];
+LightPurple=[194/255 165/255 207/255];
+LightGreen=[166/255 219/255 160/255];
+DarkGreen=[77/255 172/255 38/255];
+DarkPink=[208/255 28/255 139/255];
+MyBlue=[0.4940    0.1840    0.5560];
+MyOrange=[0.9290    0.6940    0.1250];
+DS=[51/255 160/255 44/255];
+Pav=[180/255 41/255 230/255];
+
+BSIZE=0.01;
+Dura=[-0.5 1]; Tm=Dura(1):BSIZE:Dura(2);
+Baseline=[-2 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %used to calculate Z-scores
+Tbin=-1:0.005:1; %window used to determine the optimal binsize
+PStat=0.05;
+prewin=[Dura(1) 0];
+postwin=[.1 .3];
+postwin2=[Dura(1):0.05:Dura(2)]; %bounds should match Dura
+Slopebounds=[-.5:.05:.5];
+R2Bounds=[0:0.05:0.5];
+PvalBounds=[0:0.01:0.5];
+Struct=10;
+
+%% Figure 1 VP neurons encode cue value in both tasks
+
+BrainRegion=1; %10 for Core, 100 for dms and 1000 for shell
+LeftEv=1;          RightEv=2;
+ExcInh=1;        % 1 for excitations, 2 for inhibitions -- Used to find Exc or Inh onsets in R.Ev(.).Onsets
+
+Pstat=0.05;
+stepsize=0.3;      % Used for the scatter plot projections
+
+%For COLOR-CODED PSTH only (not Average PSTH)
+Norm=0;          % 0 = no norm - max accross all the selected neurons
+% 1 each neuron is normalized wh its own max
+% 2 for baseline substraction
+
+if LeftEv==1, Xaxis=[-0.25 .5]; elseif LeftEv==3, Xaxis=[-0.25 .5]; elseif LeftEv==5; Xaxis=[-3 3]; elseif LeftEv==7; Xaxis=[-1 2]; elseif LeftEv==9; Xaxis=[-0.25 1]; end
+if ExcInh==1, Yaxis=[-5 10]; elseif ExcInh==2, Yaxis=[-5 5]; end
+
+%----------- COLORS ---------
+myblue=[0 113/255 188/255];
+mypurple=[200/255 0 255/255];
+mydarkblue=[0 0/255 255/255];
+TickSize=[0.015 0.02];
+
+DarkPurple=[123/255 50/255 148/255];
+LightPurple=[194/255 165/255 207/255];
+LightGreen=[166/255 219/255 160/255];
+DarkGreen=[77/255 172/255 38/255];
+DarkPink=[208/255 28/255 139/255];
+
+
+
+Ishow=find(RvppasALL.Param.Tm>=Xaxis(1) & RvppasALL.Param.Tm<=Xaxis(2));
+time=RvppasALL.Param.Tm(Ishow);
+if ExcInh==1, A=1; elseif ExcInh==2; A=-1;end
+figure;
+for i=1:2 %plots 3 different conditions (eg LP+No0 / LP+No+ / LP0No+)
+    
+    selectionInst=RvpppasALL.Structure==10 & ~isnan(RvpppasALL.Ev(LeftEv).Meanz) & ~isnan(RvpppasALL.Ev(RightEv).Meanz) & RvpppasALL.CSPlusRatio>=.5 & (RvpppasALL.CSPlusRatio./(RvpppasALL.CSPlusRatio+RvpppasALL.CSMinusRatio))>=.7;% & strcmp(RvpppasALL.Ninfo(:,1),'PA168Sess13_sorted-01.nex');%'PA167Sess07_sorted.nex''PA168Sess13_sorted-01.nex' & RvpppasALL.Ev(1).RespDir==0;% & CvpppasALL.STATSpost>=0.05;
+    selectionPav=RvppasALL.Structure==10 & ~isnan(RvppasALL.Ev(LeftEv).Meanz) & ~isnan(RvppasALL.Ev(RightEv).Meanz) & RvppasALL.CSPlusRatio>=.5 & (RvppasALL.CSPlusRatio./(RvppasALL.CSPlusRatio+RvppasALL.CSMinusRatio))>=.7;% & RvppasALL.Ev(1).RespDir==0;% & CvppasALL.STATSpost>=0.05;
+    
+    %------------- SORTING -------------
+    SortEv=LeftEv; NonSortEv=RightEv;
+    if i==1;
+        selection=selectionInst;
+        TMP=-RvpppasALL.Ev(SortEv).Meanz(selection,ExcInh);
+        %TMP=RvpppasALL.Ev(SortEv).RespDir(selection,ExcInh);
+    else;
+        selection=selectionPav;
+        TMP=-RvppasALL.Ev(SortEv).Meanz(selection,ExcInh);
+        %TMP=RvppasALL.Ev(SortEv).RespDir(selection,ExcInh);
+    end
+    
+    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
+    [~,SORTimg]=sort(TMP);
+    
+    %NNids are the indexes for the 3 different selections
+    if i==1; NNid1=cell2mat(RvpppasALL.Ninfo(selection,3)); NNid1=NNid1(SORTimg); end
+    if i==2; NNid2=cell2mat(RvppasALL.Ninfo(selection,3)); NNid2=NNid2(SORTimg); end
+    
+    if i==1,
+        selection=selectionInst; SortEv=LeftEv; NonSortEv=RightEv;
+        % ------------- NORMALIZATION -------------
+        if Norm==0
+            imgLeft=RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow);
+            imgRight=RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow);
+            E=[-200 400];Clim=sign(E).*abs(E).^(1/4);%ColorMap
+        elseif Norm~=0
+            imgLeft=normalize(RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow),RvpppasALL.Ev(LeftEv).PSTHz(selection,(RvpppasALL.Param.Tm<0)),Norm); %normalize to max response
+            imgRight=normalize(RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow),RvpppasALL.Ev(LeftEv).PSTHz(selection,(RvpppasALL.Param.Tm<0)),Norm); %normalize to max response of the left event
+            if Norm==1, Clim=[0 1]; elseif Norm==2, Clim=[-5 10]; end
+        end
+        
+        imgLeft=imgLeft(SORTimg,:);
+        imgRight=imgRight(SORTimg,:);
+        
+    else if i==2,
+            selection=selectionPav; SortEv=LeftEv; NonSortEv=RightEv;
+            if Norm==0
+                imgLeft=RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow);
+                imgRight=RvppasALL.Ev(RightEv).PSTHz(selection,Ishow);
+                E=[-200 400];Clim=sign(E).*abs(E).^(1/4);%ColorMap
+            elseif Norm~=0
+                imgLeft=normalize(RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow),RvppasALL.Ev(LeftEv).PSTHz(selection,(RvppasALL.Param.Tm<0)),Norm); %normalize to max response
+                imgRight=normalize(RvppasALL.Ev(RightEv).PSTHz(selection,Ishow),RvppasALL.Ev(LeftEv).PSTHz(selection,(RvppasALL.Param.Tm<0)),Norm); %normalize to max response of the left event
+                if Norm==1, Clim=[0 1]; elseif Norm==2, Clim=[-5 10]; end
+            end
+            
+            imgLeft=imgLeft(SORTimg,:);
+            imgRight=imgRight(SORTimg,:);
+            
+        end
+    end
+    
+    %---------------- HEAT MAPS ----------------------
+    subplot(3,4,i+((i-1)*3));
+    imagesc(time,[1 size(imgLeft,1)],imgLeft,Clim); %colorbar;axis tight;
+    colormap(parula);
+    ylabel('Neuron #');
+    xlabel('Time from Cue (sec)');
+    if i==1;
+        title('DS');
+    elseif i==2;
+        title('CS+');
+    end
+    
+    subplot(3,4,i+1+((i-1)*3));
+    imagesc(time,[1 size(imgRight,1)],imgRight,Clim) ;%colorbar;axis tight;
+    colormap(parula);
+    %title(Erefnames(RightEv));
+    ylabel('Neuron #');
+    xlabel('Time from Cue (sec)');
+    if i==1;
+        title('NS');
+    elseif i==2;
+        title('CS-');
+    end
+    
+    subplot(3,2,i*2);
+    hold on;
+    if i==1 %PPS
+        if     LeftEv==1, Stat=RvpppasALL.CSStat;
+        elseif LeftEv==3, Stat=RvpppasALL.CSPlusResponseStat;
+        end
+        selection1=selectionInst & Stat>=0.05;
+        selection2=selectionInst & Stat<0.05 & RvpppasALL.Ev(LeftEv).Meanz>RvpppasALL.Ev(RightEv).Meanz;
+        selection3=selectionInst & Stat<0.05 & RvpppasALL.Ev(LeftEv).Meanz<RvpppasALL.Ev(RightEv).Meanz;
+        scatter(RvpppasALL.Ev(LeftEv).Meanz(selection1),RvpppasALL.Ev(RightEv).Meanz(selection1),30, 'k');%,'filled');
+        scatter(RvpppasALL.Ev(LeftEv).Meanz(selection2),RvpppasALL.Ev(RightEv).Meanz(selection2),30, DarkPink);
+        scatter(RvpppasALL.Ev(LeftEv).Meanz(selection3),RvpppasALL.Ev(RightEv).Meanz(selection3),30, DarkGreen);
+        %selection4= R.Ev(RightEv).RespDir==1 & R.Ev(6).RespDir==1;
+        %scatter(R.Ev(LeftEv).Meanz(selection4),R.Ev(RightEv).Meanz(selection4),30,'g','filled');
+        xlabel('DS Response Magnitude (z-score)');
+        ylabel('NS Response magnitude (z-score)');
+        legend('DS=NS','DS>NS','DS<NS','Location','northwest');
+    else %PS
+        if     LeftEv==1, Stat=RvppasALL.CSStat;
+        elseif LeftEv==3, Stat=RvppasALL.CSResponseStat;
+        end
+        selection1=selectionPav & Stat>=0.05;
+        selection2=selectionPav & Stat<0.05 & RvppasALL.Ev(LeftEv).Meanz>RvppasALL.Ev(RightEv).Meanz;
+        selection3=selectionPav & Stat<0.05 & RvppasALL.Ev(LeftEv).Meanz<RvppasALL.Ev(RightEv).Meanz;
+        scatter(RvppasALL.Ev(LeftEv).Meanz(selection1),RvppasALL.Ev(RightEv).Meanz(selection1),30, 'k');%,'filled');
+        scatter(RvppasALL.Ev(LeftEv).Meanz(selection2),RvppasALL.Ev(RightEv).Meanz(selection2),30, DarkPink);
+        scatter(RvppasALL.Ev(LeftEv).Meanz(selection3),RvppasALL.Ev(RightEv).Meanz(selection3),30, DarkGreen);
+        Min=floor(min([RvppasALL.Ev(LeftEv).Meanz(selection2);RvppasALL.Ev(RightEv).Meanz(selection2)]));
+        Max=ceil(max([RvppasALL.Ev(LeftEv).Meanz(selection2);RvppasALL.Ev(RightEv).Meanz(selection2)]));
+        xlabel('CS+ Response Magnitude (z-score)');
+        ylabel('CS- Response magnitude (z-score)');
+        legend('CS+=CS-','CS+>CS-','CS+<CS-','Location','northwest');
+    end
+    Min=floor(min([RvpppasALL.Ev(LeftEv).Meanz(RvpppasALL.Structure==10);RvpppasALL.Ev(RightEv).Meanz(RvpppasALL.Structure==10)]));
+    Max=ceil(max([RvpppasALL.Ev(LeftEv).Meanz(RvpppasALL.Structure==10);RvpppasALL.Ev(RightEv).Meanz(RvpppasALL.Structure==10)]));
+    line([Min Max],[Min Max],'Color', 'k');
+    line([Min Max], [0 0],'Color', 'k');
+    line([0 0], [Min Max],'Color', 'k');
+    
+    
+    subplot(10,6,58); %Histograms for scatterplots
+    if     LeftEv==1, Stat=RvpppasALL.CSStat;
+    elseif LeftEv==3, Stat=RvpppasALL.CSPlusResponseStat;
+    end
+    axis([-10 10 0 .25]);
+    edges = [-10:1:9];
+    selection1=selectionInst & Stat>=0.05;
+    selection2=selectionInst & Stat<0.05 & RvpppasALL.Ev(LeftEv).Meanz>RvpppasALL.Ev(RightEv).Meanz;
+    selection3=selectionInst & Stat<0.05 & RvpppasALL.Ev(LeftEv).Meanz<RvpppasALL.Ev(RightEv).Meanz;
+    Delta1= RvpppasALL.Ev(LeftEv).Meanz(selection1)-RvpppasALL.Ev(RightEv).Meanz(selection1);
+    Delta2= RvpppasALL.Ev(LeftEv).Meanz(selection2)-RvpppasALL.Ev(RightEv).Meanz(selection2);
+    Delta3= RvpppasALL.Ev(LeftEv).Meanz(selection3)-RvpppasALL.Ev(RightEv).Meanz(selection3);
+    Delta1Hist=hist(Delta1,edges)/314;
+    Delta2Hist=hist(Delta2,edges)/314;
+    Delta3Hist=hist(Delta3,edges)/314;
+    edges = [-10:1:10];
+    histogram('FaceColor','k', 'BinEdges',edges,'BinCounts',Delta1Hist);
+    hold on;
+    histogram('FaceColor',DarkPink, 'BinEdges',edges,'BinCounts',Delta2Hist);
+    histogram('FaceColor',DarkGreen, 'BinEdges',edges,'BinCounts',Delta3Hist);
+    
+    
+    subplot(10,6,59); %Histograms for scatterplots
+    if     LeftEv==1, Stat=RvppasALL.CSStat;
+    elseif LeftEv==3, Stat=RvppasALL.CSrespStat;
+    end
+    axis([-10 10 0 .25]);
+    edges = [-10:1:9];
+    selection1=selectionPav & Stat>=0.05;
+    selection2=selectionPav & Stat<0.05 & RvppasALL.Ev(LeftEv).Meanz>RvppasALL.Ev(RightEv).Meanz;
+    selection3=selectionPav & Stat<0.05 & RvppasALL.Ev(LeftEv).Meanz<RvppasALL.Ev(RightEv).Meanz;
+    Delta1= RvppasALL.Ev(LeftEv).Meanz(selection1)-RvppasALL.Ev(RightEv).Meanz(selection1);
+    Delta2= RvppasALL.Ev(LeftEv).Meanz(selection2)-RvppasALL.Ev(RightEv).Meanz(selection2);
+    Delta3= RvppasALL.Ev(LeftEv).Meanz(selection3)-RvppasALL.Ev(RightEv).Meanz(selection3);
+    Delta1Hist=hist(Delta1,edges)/392;
+    Delta2Hist=hist(Delta2,edges)/392;
+    Delta3Hist=hist(Delta3,edges)/392;
+    edges = [-10:1:10];
+    histogram('FaceColor','k', 'BinEdges',edges,'BinCounts',Delta1Hist);
+    hold on;
+    histogram('FaceColor',DarkPink, 'BinEdges',edges,'BinCounts',Delta2Hist);
+    histogram('FaceColor',DarkGreen, 'BinEdges',edges,'BinCounts',Delta3Hist);
+    
+    selection1=selectionInst & ~isnan(RvpppasALL.Ev(LeftEv).Meanz) & ~isnan(RvpppasALL.Ev(RightEv).Meanz)& RvpppasALL.Ev(1).RespDir==1;% & CvpppasALL.STATSpost>=0.05;
+    selection2=selectionPav & ~isnan(RvppasALL.Ev(LeftEv).Meanz) & ~isnan(RvppasALL.Ev(RightEv).Meanz)& RvppasALL.Ev(1).RespDir==1;% & CvppasALL.STATSpost>=0.05;
+    
+    %subplot of DS vs NS excited neurons
+    if i==1
+        selection=selection1;
+        subplot(6,4,i+16);
+        psthLeft=nanmean(RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
+        semLeft=nanste(RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
+        upLeft=psthLeft+semLeft;
+        downLeft=psthLeft-semLeft;
+        hold on;
+        
+        plot(time,psthLeft,'Color',DS,'linewidth',1.5);
+        
+        psthRight=nanmean(RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
+        semRight=nanste(RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
+        upRight=psthRight+semRight;
+        downRight=psthRight-semRight;
+        
+        plot(time,psthRight,'k','linewidth',1.5);
+        %legend('DS','NS','Location','northwest');
+        patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],DS,'EdgeColor','none');alpha(0.5);
+        patch([time,time(end:-1:1)],[upRight,downRight(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
+        plot([0 0], [-10 20],'LineStyle',':','Color','k');
+        plot([-10 20], [0 0],'LineStyle',':','Color','k');
+        if LeftEv==3
+            plot([1 1], [-10 20],'LineStyle',':','Color','k');
+            plot([2.5 2.5], [-10 20],'LineStyle',':','Color','k');
+        elseif LeftEv==5
+            plot([-2 -2], [-10 20],'LineStyle',':','Color','k');
+            plot([.5 .5], [-10 20],'LineStyle',':','Color','k');
+        elseif LeftEv==7
+            plot([1 1], [-10 20],'LineStyle',':','Color','k');
+            plot([5 5], [-10 20],'LineStyle',':','Color','k');
+        end
+        axis([Xaxis,Yaxis]);
+        %xlabel('Time from Cue (sec)');
+        ylabel('Average Response (z-score)');
+        
+    else if i==2
+            selection=selection2;
+            subplot(6,4,i+19);
+            psthLeft=nanmean(RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
+            semLeft=nanste(RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
+            upLeft=psthLeft+semLeft;
+            downLeft=psthLeft-semLeft;
+            hold on;
+            
+            plot(time,psthLeft,'Color',Pav,'linewidth',1.5);
+            
+            psthRight=nanmean(RvppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
+            semRight=nanste(RvppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
+            upRight=psthRight+semRight;
+            downRight=psthRight-semRight;
+            
+            plot(time,psthRight,'k','linewidth',1.5);
+            %legend('CS+','CS-','Location','northwest');
+            patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],Pav,'EdgeColor','none');alpha(0.5);
+            patch([time,time(end:-1:1)],[upRight,downRight(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
+            plot([0 0], [-10 20],'LineStyle',':','Color','k');
+            plot([-10 20], [0 0],'LineStyle',':','Color','k');
+            if LeftEv==3
+                plot([1 1], [-10 20],'LineStyle',':','Color','k');
+                plot([2.5 2.5], [-10 20],'LineStyle',':','Color','k');
+            elseif LeftEv==5
+                plot([-2 -2], [-10 20],'LineStyle',':','Color','k');
+                plot([.5 .5], [-10 20],'LineStyle',':','Color','k');
+            elseif LeftEv==7
+                plot([1 1], [-10 20],'LineStyle',':','Color','k');
+                plot([5 5], [-10 20],'LineStyle',':','Color','k');
+            end
+            axis([Xaxis,Yaxis]);
+            
+            xlabel('Time from Cue (sec)');
+            ylabel('Average Response (z-score)');
+        end
+    end
+        
+        
+        selection1=selectionInst & ~isnan(RvpppasALL.Ev(LeftEv).Meanz) & ~isnan(RvpppasALL.Ev(RightEv).Meanz)& RvpppasALL.Ev(1).RespDir==-1;% & CvpppasALL.STATSpost>=0.05;
+        selection2=selectionPav & ~isnan(RvppasALL.Ev(LeftEv).Meanz) & ~isnan(RvppasALL.Ev(RightEv).Meanz)& RvppasALL.Ev(1).RespDir==-1;% & CvppasALL.STATSpost>=0.05;
+        
+        if i==1
+            selection=selection1;
+            subplot(6,4,i+17);
+            psthLeft=nanmean(RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
+            semLeft=nanste(RvpppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
+            upLeft=psthLeft+semLeft;
+            downLeft=psthLeft-semLeft;
+            hold on;
+            
+            plot(time,psthLeft,'Color',DS,'linewidth',1.5);
+            
+            psthRight=nanmean(RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
+            semRight=nanste(RvpppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
+            upRight=psthRight+semRight;
+            downRight=psthRight-semRight;
+            
+            plot(time,psthRight,'k','linewidth',1.5);
+            legend('DS','NS','Location','northwest');
+            patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],DS,'EdgeColor','none');alpha(0.5);
+            patch([time,time(end:-1:1)],[upRight,downRight(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
+            plot([0 0], [-10 20],'LineStyle',':','Color','k');
+            plot([-10 20], [0 0],'LineStyle',':','Color','k');
+            if LeftEv==3
+                plot([1 1], [-10 20],'LineStyle',':','Color','k');
+                plot([2.5 2.5], [-10 20],'LineStyle',':','Color','k');
+            elseif LeftEv==5
+                plot([-2 -2], [-10 20],'LineStyle',':','Color','k');
+                plot([.5 .5], [-10 20],'LineStyle',':','Color','k');
+            elseif LeftEv==7
+                plot([1 1], [-10 20],'LineStyle',':','Color','k');
+                plot([5 5], [-10 20],'LineStyle',':','Color','k');
+            end
+            axis([Xaxis,Yaxis]);
+            
+            %xlabel('Time from Cue (sec)');
+            %ylabel('Average Response (z-score)');
+            
+        else if i==2
+                selection=selection2;
+                subplot(6,4,i+20);
+                psthLeft=nanmean(RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
+                semLeft=nanste(RvppasALL.Ev(LeftEv).PSTHz(selection,Ishow),1);
+                upLeft=psthLeft+semLeft;
+                downLeft=psthLeft-semLeft;
+                hold on;
+                
+                plot(time,psthLeft,'Color',Pav,'linewidth',1.5);
+                
+                psthRight=nanmean(RvppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
+                semRight=nanste(RvppasALL.Ev(RightEv).PSTHz(selection,Ishow),1);
+                upRight=psthRight+semRight;
+                downRight=psthRight-semRight;
+                
+                plot(time,psthRight,'k','linewidth',1.5);
+                legend('CS+','CS-','Location','northwest');
+                patch([time,time(end:-1:1)],[upLeft,downLeft(end:-1:1)],Pav,'EdgeColor','none');alpha(0.5);
+                patch([time,time(end:-1:1)],[upRight,downRight(end:-1:1)],'k','EdgeColor','none');alpha(0.5);
+                plot([0 0], [-10 20],'LineStyle',':','Color','k');
+                plot([-10 20], [0 0],'LineStyle',':','Color','k');
+                if LeftEv==3
+                    plot([1 1], [-10 20],'LineStyle',':','Color','k');
+                    plot([2.5 2.5], [-10 20],'LineStyle',':','Color','k');
+                elseif LeftEv==5
+                    plot([-2 -2], [-10 20],'LineStyle',':','Color','k');
+                    plot([.5 .5], [-10 20],'LineStyle',':','Color','k');
+                elseif LeftEv==7
+                    plot([1 1], [-10 20],'LineStyle',':','Color','k');
+                    plot([5 5], [-10 20],'LineStyle',':','Color','k');
+                end
+                axis([Xaxis,Yaxis]);
+                
+                xlabel('Time from Cue (sec)');
+                %ylabel('Average Response (z-score)');
+            end
+        end
+    
+end
+
+subplot(9,2,[14 16]);
+%Plot of auROC descriptions for CS versus DS on response versus no response
+%(percent of neurons!)
+sel1=RvpppasALL.Structure~=0;
+sel2=RvppasALL.Structure~=0;
+axis([0 1 0 .25]);
+%histogram(-RvpppasALL.CSauROC(RvpppasALL.Structure==10),[0:0.05:1],'Normalization','probability');
+cdfplot(-RvpppasALL.CSauROC(RvpppasALL.Structure==10 & RvpppasALL.CSPlusRatio>=.5 & (RvpppasALL.CSPlusRatio./(RvpppasALL.CSPlusRatio+RvpppasALL.CSMinusRatio))>=.7));
+hold on;
+cdfplot(-RvppasALL.CSauROC(RvppasALL.Structure==10 & RvppasALL.CSPlusRatio>=.5 & (RvppasALL.CSPlusRatio./(RvppasALL.CSPlusRatio+RvppasALL.CSMinusRatio))>=.7));
+%histogram(-RvppasALL.CSauROC(RvppasALL.Structure==10),[0:0.05:1],'Normalization','probability');
+MEANsel1=nanmedian(-RvpppasALL.CSauROC(RvpppasALL.Structure==10 & RvpppasALL.CSPlusRatio>=.5 & (RvpppasALL.CSPlusRatio./(RvpppasALL.CSPlusRatio+RvpppasALL.CSMinusRatio))>=.7));
+MEANsel2=nanmedian(-RvppasALL.CSauROC(RvppasALL.Structure==10 & RvppasALL.CSPlusRatio>=.5 & (RvppasALL.CSPlusRatio./(RvppasALL.CSPlusRatio+RvppasALL.CSMinusRatio))>=.7));
+plot([MEANsel1 MEANsel1], [0 1],'Color',DS);
+plot([MEANsel2 MEANsel2], [0 1],'Color',Pav);
+plot([.5 .5], [0 1],'Color','k');
+legend('DS','CS+','Location','northwest');
+xlabel('Cue Identity auROC');
+ylabel('Proportion of Neurons');
+
+task=cat(1,ones(length(RvpppasALL.CSauROC(selectionInst)),1),zeros(length(RvppasALL.CSauROC(selectionPav)),1));
+rat=cat(1,RvpppasALL.Rat(selectionInst),RvppasALL.Rat(selectionPav));
+tbl=table(cat(1,-RvpppasALL.CSauROC(selectionInst),-RvppasALL.CSauROC(selectionPav)),task,rat,'variablenames',{'auROC','Task','Subject'});
+%lme=fitlme(tbl,'auROC~Task+(1|Subject)');%+(1|Subject:Exposure)');
+
+
+%stats on firing of excited and inhibited neurons
+excitedNsP=selectionPav & (RvppasALL.Ev(1).RFLAG==1 | RvppasALL.Ev(1).RespDir==1) & RvppasALL.Ev(1).RFLAG~=-1 & RvppasALL.Ev(1).RespDir~=-1;
+excitedNsI=selectionInst & (RvpppasALL.Ev(1).RFLAG==1 | RvpppasALL.Ev(1).RespDir==1) & RvpppasALL.Ev(1).RFLAG~=-1 & RvpppasALL.Ev(1).RespDir~=-1;
+inhibNsI=selectionInst & (RvpppasALL.Ev(1).RFLAG==-1 | RvpppasALL.Ev(1).RespDir==-1) & RvpppasALL.Ev(1).RFLAG~=1 & RvpppasALL.Ev(1).RespDir~=1;
+inhibNsP=selectionPav & (RvppasALL.Ev(1).RFLAG==-1 | RvppasALL.Ev(1).RespDir==-1) & RvppasALL.Ev(1).RFLAG~=1 & RvppasALL.Ev(1).RespDir~=1;
+
+% [~,p,~,stats]=ttest(RvppasALL.Ev(1).Meanraw(excitedNsP),RvppasALL.Ev(2).Meanraw(excitedNsP)); %CS excited
+% [~,p,~,stats]=ttest(RvppasALL.Ev(1).Meanraw(inhibNsP),RvppasALL.Ev(2).Meanraw(inhibNsP)); %CS inhibited
+% [~,p,~,stats]=ttest(RvpppasALL.Ev(1).Meanraw(excitedNsI),RvpppasALL.Ev(2).Meanraw(excitedNsI)); %DS excited
+% [~,p,~,stats]=ttest(RvpppasALL.Ev(1).Meanraw(inhibNsI),RvpppasALL.Ev(2).Meanraw(inhibNsI)); %DS inhibited
+%---------------- COLOR SCALE ----------------------
+subplot(10,6,60);
+imagesc(time,[1 size(imgLeft,1)],imgLeft,Clim);colorbar;
+%ylabel('z-score');
+
+%
+% if SAVEFIG==1
+%     Figpath='C:\Users\jricha95\Documents\Documents (3)\Experiments\New Experiments\PS PPS';
+%     exportfig(figure(1), sprintf('%s\\Figure1cmyk.eps',Figpath), 'FontMode', 'fixed','FontSize', 12,'LineWidth', 0.5, 'Width', 10, 'format','eps2','color','cmyk','renderer', 'painters','resolution',300);
+% %     %end
+
+
+
+
+
+%% Figure 5 pie charts
+
+figure;
+
+
+%Pie charts of response patterns
+%DS
+subplot(4,5,11);
+exc=sum((RvpppasALL.Ev(1).RFLAG(selectionInst)==1 | RvpppasALL.Ev(1).RespDir(selectionInst)==1) & RvpppasALL.Ev(1).RFLAG(selectionInst)~=-1 & RvpppasALL.Ev(1).RespDir(selectionInst)~=-1);
+inh=sum(RvpppasALL.Ev(1).RespDir(selectionInst)==-1 | RvpppasALL.Ev(1).RFLAG(selectionInst)==-1);
+no=sum(RvpppasALL.Ev(1).RFLAG(selectionInst)==0 & RvpppasALL.Ev(1).RespDir(selectionInst)==0);
+response=[exc inh no];
+labels = {'Excited','Inhibited','No Response'};
+pie(response);
+colormap([252/255 206/255 46/255;      %// excitation yellow
+    50/255 67/255 186/255;      %// inhibition blue
+    46/255 183/255 164/255]);      %// no response teal
+title('DS');
+
+%CS+
+subplot(4,5,12);
+exc=sum((RvppasALL.Ev(1).RFLAG(selectionPav)==1 | RvppasALL.Ev(1).RespDir(selectionPav)==1) & RvppasALL.Ev(1).RFLAG(selectionPav)~=-1 & RvppasALL.Ev(1).RespDir(selectionPav)~=-1);
+inh=sum(RvppasALL.Ev(1).RespDir(selectionPav)==-1 | RvppasALL.Ev(1).RFLAG(selectionPav)==-1);
+no=sum(RvppasALL.Ev(1).RFLAG(selectionPav)==0 & RvppasALL.Ev(1).RespDir(selectionPav)==0);
+response=[exc inh no];
+pie(response);
+title('CS+');
+
+%DS port entry
+subplot(4,5,13);
+exc=sum((RvpppasALL.Ev(8).RFLAG(selectionInst)==1 | RvpppasALL.Ev(8).RespDir(selectionInst)==1) & RvpppasALL.Ev(5).RFLAG(selectionInst)~=-1 & RvpppasALL.Ev(5).RespDir(selectionInst)~=-1);
+inh=sum(RvpppasALL.Ev(8).RespDir(selectionInst)==-1 | RvpppasALL.Ev(8).RFLAG(selectionInst)==-1);
+no=sum(RvpppasALL.Ev(8).RFLAG(selectionInst)==0 & RvpppasALL.Ev(8).RespDir(selectionInst)==0);
+response=[exc inh no];
+labels = {'Excited','Inhibited','No Response'};
+pie(response);
+title('post-DS Port Entry');
+
+%CS+ port entry
+subplot(4,5,14);
+exc=sum((RvppasALL.Ev(8).RFLAG(selectionPav)==1 | RvppasALL.Ev(8).RespDir(selectionPav)==1) & RvppasALL.Ev(5).RFLAG(selectionPav)~=-1 & RvppasALL.Ev(5).RespDir(selectionPav)~=-1);
+inh=sum(RvppasALL.Ev(8).RespDir(selectionPav)==-1 | RvppasALL.Ev(8).RFLAG(selectionPav)==-1);
+no=sum(RvppasALL.Ev(8).RFLAG(selectionPav)==0 & RvppasALL.Ev(8).RespDir(selectionPav)==0);
+response=[exc inh no];
+h=pie(response);
+title('post-CS+ Port Entry');
+
+%CS+ reward
+subplot(4,5,15);
+exc=sum((RvppasALL.Ev(9).RFLAG(selectionPav)==1 | RvppasALL.Ev(9).RespDir(selectionPav)==1) & RvppasALL.Ev(7).RFLAG(selectionPav)~=-1 & RvppasALL.Ev(7).RespDir(selectionPav)~=-1);
+inh=sum(RvppasALL.Ev(9).RespDir(selectionPav)==-1 | RvppasALL.Ev(9).RFLAG(selectionPav)==-1);
+no=sum(RvppasALL.Ev(9).RFLAG(selectionPav)==0 & RvppasALL.Ev(9).RespDir(selectionPav)==0);
+response=[exc inh no];
+h=pie(response);
+title('post-CS+ Reward');
+
+%% Figure 5 behavior
+clear all;
+
+load('RAWvppsALL.mat')
+load('RAWvpppsALL.mat')
+load('RAWvppasALL.mat')
+load('RAWvpppasALL.mat')
+load('RPPSall.mat')
+load('RPSall.mat')
+load('RvppasALL.mat')
+load('RvpppasALL.mat')
+
+
+SessRat=[];
+for task=1:4
+    if task==1 %inst suc
+        RAW=RAWvpppsALL;
+        R=RPPSall;
+    elseif task==2 %inst suc(alc)
+        RAW=RAWvpppasALL;
+        R=RvpppasALL;
+    elseif task==3 %pav suc
+        RAW=RAWvppsALL;
+        R=RPSall;
+    elseif task==4 %pav suc(alc)
+        RAW=RAWvppasALL;
+        R=RvppasALL;
+    end
+    
+    %which neurons are included
+    selection=R.Structure==10 & R.CSPlusRatio>=.5 & (R.CSPlusRatio./(R.CSPlusRatio+R.CSMinusRatio))>=.7;
+    
+    %which sessions are included
+    allsessions=unique(R.Ninfo(:,1)); %all sessions from this task
+    [includedsessions,index]=unique(R.Ninfo(selection,1)); %all included sessions
+    included=ismember(allsessions,includedsessions); %which of all sessions are included
+    
+    %get rat for each session
+    Rats=R.Rat(index)+task*100; %make sure rats with same name from different experiment aren't accidentally combined
+    SessRat=cat(1,SessRat,Rats,Rats); %twice to represent the CSPlus data and then the CSMinus data
+
+    
+    
+    RAW=RAW(included);
+    
+    for i=1:length(RAW)
+        if task==1 | task==3
+            Ev1=strcmp('CSPlus1', RAW(i).Einfo(:,2));
+            Ev2=strcmp('CSMinus1', RAW(i).Einfo(:,2));
+        else
+            Ev1=strcmp('CSPlus', RAW(i).Einfo(:,2));
+            Ev2=strcmp('CSMinus', RAW(i).Einfo(:,2));
+        end
+        PE=strcmp('PortEntry', RAW(i).Einfo(:,2));
+        
+        CSPluses=RAW(i).Erast{Ev1};
+        CSMinuses=RAW(i).Erast{Ev2};
+        PortEntries=RAW(i).Erast{PE};
+        
+        CSPLatency=[];
+        for j=1:length(CSPluses)
+            if sum(PortEntries>CSPluses(j))>0
+                CSPLatency(j,1)=min(PortEntries(PortEntries>CSPluses(j)))-CSPluses(j);
+            else
+                CSPLatency(j,1)=100;
+            end
+        end
+        CSPLatency(CSPLatency>10)=[];
+        
+        CSMLatency=[];
+        for j=1:length(CSMinuses)
+            if sum(PortEntries>CSMinuses(j))>0
+                CSMLatency(j,1)=min(PortEntries(PortEntries>CSMinuses(j)))-CSMinuses(j);
+            else
+                CSMLatency(j,1)=100; %will be removed
+            end
+        end
+        CSMLatency(CSMLatency>10)=[];        
+        
+        CSPMean{task,1}(i,1)=mean(CSPLatency,'omitnan');
+        CSMMean{task,1}(i,1)=mean(CSMLatency,'omitnan');
+        CSPlusProb{task,1}(i,1)=length(CSPLatency)/length(CSPluses);
+        CSMinusProb{task,1}(i,1)=length(CSMLatency)/length(CSMinuses);
+        
+    end
+    
+end
+
+%% plotting
+
+figure;
+
+DS=[0    0.4470    0.7410];
+Pav=[0.8500    0.3250    0.0980];
+DSAlcohol=[51/255 160/255 44/255];
+PavAlcohol=[180/255 41/255 230/255];
+black=[0 0 0];
+Colors=[DS;black;DSAlcohol;black;Pav;black;PavAlcohol;black];
+
+
+%probability
+data=[CSPlusProb{1,1};CSMinusProb{1,1};CSPlusProb{2,1};CSMinusProb{2,1};CSPlusProb{3,1};CSMinusProb{3,1};CSPlusProb{4,1};CSMinusProb{4,1}];
+group=[ones(size(CSPlusProb{1,1}));2*ones(size(CSMinusProb{1,1}));3*ones(size(CSPlusProb{2,1}));4*ones(size(CSMinusProb{2,1}));5*ones(size(CSPlusProb{3,1}));6*ones(size(CSMinusProb{3,1}));7*ones(size(CSPlusProb{4,1}));8*ones(size(CSMinusProb{4,1}))];
+task=[ones(size(CSPlusProb{1,1}));ones(size(CSMinusProb{1,1}));ones(size(CSPlusProb{2,1}));ones(size(CSMinusProb{2,1}));2*ones(size(CSPlusProb{3,1}));2*ones(size(CSMinusProb{3,1}));2*ones(size(CSPlusProb{4,1}));2*ones(size(CSMinusProb{4,1}))];
+exposure=[ones(size(CSPlusProb{1,1}));ones(size(CSMinusProb{1,1}));2*ones(size(CSPlusProb{2,1}));2*ones(size(CSMinusProb{2,1}));ones(size(CSPlusProb{3,1}));ones(size(CSMinusProb{3,1}));2*ones(size(CSPlusProb{4,1}));2*ones(size(CSMinusProb{4,1}))];
+cue=[ones(size(CSPlusProb{1,1}));2*ones(size(CSMinusProb{1,1}));ones(size(CSPlusProb{2,1}));2*ones(size(CSMinusProb{2,1}));ones(size(CSPlusProb{3,1}));2*ones(size(CSMinusProb{3,1}));ones(size(CSPlusProb{4,1}));2*ones(size(CSMinusProb{4,1}))];
+
+subplot(1,2,1);
+hold on;
+boxplot(data,group,'colors',Colors,'symbol','');
+
+for i=1:8
+    scatter(i+((rand(sum(group==i),1)-0.5)/5),data(group==i),[],Colors(i,:),'filled')
+end
+
+plot([4.5 4.5],[0 0.9],'--','color','k');
+title('Port Entry Probability');
+xticks(1:8);
+xticklabels({'DS','NS','DS','NS','CS+','CS-','CS+','CS-'});
+text(0.8,0.6,'Sucrose','color',DS);
+text(2.8,0.8,'Suc(alc)','color',DSAlcohol);
+text(4.8,0.5,'Sucrose','color',Pav);
+text(6.8,0.6,'Suc(alc)','color',PavAlcohol);
+text(4.3,0.97,'ns','color','k');
+text(2,0.9,'<- ns ->','color','k','horizontalalignment','center');
+text(6,0.85,'<- ns ->','color','k','horizontalalignment','center');
+axis([0 9 0 1]);
+
+%stats
+%[~,ProbAnova,ProbStats]=anovan(data,{task,cue,SessRat},'varnames',{'Task','Cue','Subject'},'random',3,'nested',[0 0 0;0 0 0;1 0 0],'model','full');
+[~,ProbAnova,ProbStats]=anovan(data,{task,exposure,cue},'varnames',{'Task','Exposure','Cue'},'model','full','display','off');
+%[~,ProbAnova,ProbStats]=anovan(data(cue==1),{task(cue==1),exposure(cue==1)},'varnames',{'Task','Exposure'},'model','full','display','off');
+%multcompare(ProbStats,'dimension',[1 2]);
+
+%latency
+data=[CSPMean{1,1};CSMMean{1,1};CSPMean{2,1};CSMMean{2,1};CSPMean{3,1};CSMMean{3,1};CSPMean{4,1};CSMMean{4,1}];
+group=[ones(size(CSPMean{1,1}));2*ones(size(CSMMean{1,1}));3*ones(size(CSPMean{2,1}));4*ones(size(CSMMean{2,1}));5*ones(size(CSPMean{3,1}));6*ones(size(CSMMean{3,1}));7*ones(size(CSPMean{4,1}));8*ones(size(CSMMean{4,1}))];
+
+subplot(1,2,2);
+hold on;
+boxplot(data,group,'colors',Colors,'symbol','');
+
+for i=1:8
+    scatter(i+((rand(sum(group==i),1)-0.5)/5),data(group==i),[],Colors(i,:),'filled')
+end
+
+plot([4.5 4.5],[0 8],'--','color','k');
+title('Port Entry Latency (sec)');
+xticks(1:8);
+xticklabels({'DS','NS','DS','NS','CS+','CS-','CS+','CS-'});
+xticklabels({'DS','NS','DS','NS','CS+','CS-','CS+','CS-'});
+text(0.8,7,'Sucrose','color',DS);
+text(2.8,6,'Suc(alc)','color',DSAlcohol);
+text(4.8,7.5,'Sucrose','color',Pav);
+text(6.8,7,'Suc(alc)','color',PavAlcohol);
+text(2,1.5,'<- ns ->','color','k','horizontalalignment','center');
+text(6,1.5,'<- ns ->','color','k','horizontalalignment','center');
+text(4.3,8.5,'*','color','k');
+
+axis([0 9 0 9]);
+
+%stats
+%[~,LatAnova,LatStats]=anovan(data,{task,cue,SessRat},'varnames',{'Task','Cue','Subject'},'random',3,'nested',[0 0 0;0 0 0;1 0 0],'model','full');
+[~,LatAnova,LatStats]=anovan(data,{task,exposure,cue},'varnames',{'Task','Exposure','Cue'},'model','full','display','off');
+%[~,LatAnova,LatStats]=anovan(data(cue==1),{task(cue==1),exposure(cue==1)},'varnames',{'Task','Exposure'},'model','full','display','off');
+
+%multcompare(LatStats,'dimension',[1 2]);
+
+
+
+
+
+

+ 161 - 0
AnalysisAndFigures/m_SingleUnitDecodingAlcoholExposed.m

@@ -0,0 +1,161 @@
+%decodes trial identity from spiking across bins for each individual neuron
+clear all
+load('RAWvppsALL.mat')
+load('RAWvpppsALL.mat')
+load('RAWvppasALL.mat')
+load('RAWvpppasALL.mat')
+
+for PavInst=1:4; %DS is 1, Pav is 2
+    if PavInst==1
+        RAW=RAWvpppsALL;
+    else if PavInst==2
+            RAW=RAWvppsALL;
+        else if PavInst==3
+                RAW=RAWvpppasALL;
+            else if PavInst==4
+                    RAW=RAWvppasALL;
+                end
+            end
+        end
+    end
+    folds = 5; %number of times cross-validated:
+    Dura=[0 0.3]; %size of bin used for decoding:
+    bins = 13; %number of bins
+    binint = 0.3; %spacing of bins
+    binstart = -.9; %start time of first bin relative to event
+    shuffs = 1; %number of shuffled models created
+    trials = 30; %max number of trials to use
+    
+    
+    NN = 0;
+    
+    
+    %events being compared
+    
+    
+    UnitMdlAcc=NaN(100,bins);
+    UnitMdlAccShuff=NaN(100,bins);
+    
+    for i=1:length(RAW) %loops through sessions
+        if (PavInst==1 | PavInst==2)
+            RD1=strcmp('CSPlus1', RAW(i).Einfo(:,2));
+            RD2=strcmp('CSMinus1', RAW(i).Einfo(:,2));
+        else
+            RD1=strcmp('CSPlus', RAW(i).Einfo(:,2));
+            RD2=strcmp('CSMinus', RAW(i).Einfo(:,2));
+        end
+        
+        for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
+            NN=NN+1; %neuron counter
+            for l=1:bins
+                
+                LowHz=zeros(1,1);
+                DecodeSpikes=NaN(1,1);
+                DecodeRs=NaN(1,1);
+                
+                %find all spikes surrounding ev1
+                [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{RD1},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+                
+                %only do at most the number of trials listed above
+                if N1>trials
+                    trials1=trials;
+                else
+                    trials1=N1;
+                end
+                
+                %get the number of spikes on each trial
+                for m=1:trials1
+                    DecodeSpikes(m,1)=sum(PSR1{1,m}>(binstart));
+                    DecodeRs(m,1)=1;
+                end
+               
+                %find all spikes surrounding ev2
+                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{RD2},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+
+                %only do at most the number of trials listed above
+                if N2>trials
+                    trials2=trials;
+                else
+                    trials2=N2;
+                end
+                
+                for n=1:trials2
+                    DecodeSpikes(n+m,1)=sum(PSR2{1,n}>(binstart));
+                    DecodeRs(n+m,1)=2;
+                end
+                
+                if sum(DecodeSpikes(1:trials1,1))<7 || sum(DecodeSpikes((1+trials1):(trials1+trials2),1))<7 %prevents errors with LDA on conditions with no variance
+                    LowHz(1,1)=1;
+                end
+                
+                
+                %if one neuron or more has too few spikes, so LDA has error because not
+                %enough variance, this removes those neurons from the analysis
+                if sum(LowHz)==0
+                    
+                    %creating models
+                    CVacc = NaN(folds,1);
+                    CVaccSh = NaN(folds,1);
+                    %normal model
+                    Partitions = cvpartition(DecodeRs,'KFold',folds);
+                    for r = 1:folds
+                        
+                        LDAModel = fitcdiscr(DecodeSpikes(Partitions.training(r),:),DecodeRs(Partitions.training(r)));
+                        prediction = predict(LDAModel,DecodeSpikes(Partitions.test(r),:));
+                        actual = DecodeRs(Partitions.test(r));
+                        correct = prediction - actual;
+                        CVacc(r) = sum(correct==0) / length(correct);
+                        
+                    end
+                    
+                    %shuffled model
+                    for q=1:shuffs
+                        DecodeRsSh=DecodeRs(randperm(length(DecodeRs)));
+                        PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
+                        for s = 1:folds
+                            LDAModelSh = fitcdiscr(DecodeSpikes(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)));
+                            predictionSh = predict(LDAModelSh,DecodeSpikes(PartitionsSh.test(s),:));
+                            actualSh = DecodeRsSh(PartitionsSh.test(s));
+                            correctSh = predictionSh - actualSh;
+                            CVaccSh(s) = sum(correctSh==0) / length(correctSh);
+                        end
+                        AccShuff(q,1) = nanmean(CVaccSh);
+                    end
+                    UnitMdlAcc(NN,l) = nanmean(CVacc);
+                    UnitMdlAccShuff(NN,l) = nanmean(AccShuff);
+                else
+                    UnitMdlAcc(NN,l) = NaN;
+                    UnitMdlAccShuff(NN,l) = NaN;
+                end
+            end
+            
+            fprintf('Neuron ID # %d\n',NN);
+        end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
+    end %sessions: FOR i=1:length(RAW)
+    
+    if PavInst==1
+        UnitMdlAccDSSucrose=UnitMdlAcc;
+        UnitMdlAccDSShuffSucrose=UnitMdlAccShuff;
+        save('UnitMdlDSSucrose.mat','UnitMdlAccDSSucrose')
+        save('UnitMdlDSShuffSucrose.mat','UnitMdlAccDSShuffSucrose')
+    else if PavInst==2
+            UnitMdlAccPavSucrose=UnitMdlAcc;
+            UnitMdlAccPavShuffSucrose=UnitMdlAccShuff;
+            save('UnitMdlPavSucrose.mat','UnitMdlAccPavSucrose')
+            save('UnitMdlPavShuffSucrose.mat','UnitMdlAccPavShuffSucrose')
+        else if PavInst==3
+                UnitMdlAccDSAlcohol=UnitMdlAcc;
+                UnitMdlAccDSShuffAlcohol=UnitMdlAccShuff;
+                save('UnitMdlDSAlcohol.mat','UnitMdlAccDSAlcohol')
+                save('UnitMdlDSShuffAlcohol.mat','UnitMdlAccDSShuffAlcohol')                
+            else if PavInst==4
+                    UnitMdlAccPavAlcohol=UnitMdlAcc;
+                    UnitMdlAccPavShuffAlcohol=UnitMdlAccShuff;
+                    save('UnitMdlPavAlcohol.mat','UnitMdlAccPavAlcohol')
+                    save('UnitMdlPavShuffAlcohol.mat','UnitMdlAccPavShuffAlcohol')
+                    
+                end
+            end
+        end
+    end
+end

+ 251 - 0
AnalysisAndFigures/n_PooledDecodeAlcoholExposed.m

@@ -0,0 +1,251 @@
+%decodes trial identity from spiking across bins for each individual neuron
+
+clear all
+load('RAWvppsALL.mat')
+load('RAWvpppsALL.mat')
+load('RAWvppasALL.mat')
+load('RAWvpppasALL.mat')
+if ~exist('RPPSall'),load('RPPSall.mat'); end
+if ~exist('RPSall'),load('RPSall.mat'); end
+if ~exist('RvppasALL'),load('RvppasALL.mat'); end
+if ~exist('RvpppasALL'),load('RvpppasALL.mat'); end
+
+for PavInst=1:4; %DS is 1, Pav is 2
+    if PavInst==1
+        RAW=RAWvpppsALL;
+        selection=RPPSall.Structure==10 & RPPSall.CSPlusRatio>=.5 & (RPPSall.CSPlusRatio./(RPPSall.CSPlusRatio+RPPSall.CSMinusRatio))>=.7;
+    else if PavInst==2
+            RAW=RAWvppsALL;
+            selection=RPSall.Structure==10 & RPSall.CSPlusRatio>=.5 & (RPSall.CSPlusRatio./(RPSall.CSPlusRatio+RPSall.CSMinusRatio))>=.7;
+        else if PavInst==3
+                RAW=RAWvpppasALL;
+                selection=RvpppasALL.Structure==10 & RvpppasALL.CSPlusRatio>=.5 & (RvpppasALL.CSPlusRatio./(RvpppasALL.CSPlusRatio+RvpppasALL.CSMinusRatio))>=.7;
+            else if PavInst==4
+                    RAW=RAWvppasALL;
+                    selection=RvppasALL.Structure==10 & RvppasALL.CSPlusRatio>=.5 & (RvppasALL.CSPlusRatio./(RvppasALL.CSPlusRatio+RvppasALL.CSMinusRatio))>=.7;
+                end
+            end
+        end
+    end   
+    
+    
+    
+    
+    
+    TotalNeurons = 0;
+    
+    
+    NumNeurons=[1 5 10 25 50 100]; %matrix of how many neurons to use on each iteration
+    repetitions=20; %how many times to run the analysis
+    folds = 5; %number of times cross-validated:
+    Dura=[0 0.3]; %size of bin used for decoding:
+    bins = 1; %number of bins
+    binint = 0.3; %spacing of bins
+    binstart = 0; %start time of first bin relative to event
+    shuffs = 1; %number of shuffled models created
+    
+    
+    
+    for i=1:length(NumNeurons)
+        PoolMdlAcc{i,1}=NaN(folds,bins);
+        PoolMdlAccShuff{i,1}=NaN(folds,bins);
+    end
+    
+    %total number of neurons in all sessions and events per session
+    %for i=1:length(RAW)
+    TotalNeurons=sum(selection);
+    %     Ev1perSess(i)=length(RAW(i).Erast{Ev1,1});
+    %     Ev2perSess(i)=length(RAW(i).Erast{Ev2,1});
+    %end
+    
+    %figure out how many trials of each event can be used
+    % Ev1s=min(Ev1perSess);
+    % Ev2s=min(Ev2perSess);
+    
+    %Need to keep number of trials used in VP and NAc constant, so have to pick
+    %minimum number across all sessions in both regions. Right now with
+    %existing sessions it's 20 Ev1s and 20 Ev2s
+    Ev1s=27;
+    Ev2s=27;
+    
+    for v = 1:length(NumNeurons)
+        for u=1:repetitions
+            %pick which neurons to use
+            SetupSel=cat(1,ones(NumNeurons(v),1),zeros(TotalNeurons-NumNeurons(v),1));
+            Sel=(SetupSel(randperm(length(SetupSel)))==1);
+            
+            %setup the event identity matrix for decoding
+            DecodeRs(1:Ev1s,1)=1;
+            DecodeRs((Ev1s+1):(Ev1s+Ev2s),1)=2;
+            
+            for l=1:bins
+                DecodeSpikes=NaN(1,1);
+                AllSpikes=NaN(1,1);
+                LowHz=zeros(1,NumNeurons(v));
+                NN = 1;
+                AllNN=1;
+                for i=1:length(RAW) %loops through sessions
+                    for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
+                        if selection(AllNN)==1
+                            Ev1Spikes=NaN(1,1);
+                            Ev2Spikes=NaN(1,1);
+                            
+                            %events being compared
+                            if (PavInst==1 | PavInst==2)
+                                Ev1=strcmp('CSPlus1', RAW(i).Einfo(:,2));
+                                Ev2=strcmp('CSMinus1', RAW(i).Einfo(:,2));
+                            else
+                                Ev1=strcmp('CSPlus', RAW(i).Einfo(:,2));
+                                Ev2=strcmp('CSMinus', RAW(i).Einfo(:,2));
+                            end
+                            
+                            
+                            %get number of spikes for this neuron in this bin for all
+                            %Ev1 trials
+                            [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Ev1},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+                            for m=1:length(PSR1)
+                                Ev1Spikes(m,1)=sum(PSR1{1,m}>(binstart));
+                            end
+                            
+                            %pick which trials get used for decoding
+                            SetupTrials=cat(1,ones(Ev1s,1),zeros(length(PSR1)-Ev1s,1));
+                            Trials=SetupTrials==1; %(SetupTrials(randperm(length(SetupTrials)))==1); %randomize trials (or not)
+                            
+                            %put spikes from those trials in a matrix
+                            AllSpikes(Trials,NN)=Ev1Spikes(Trials,1);
+                            
+                            %get all the spikes from reward 1p2 trials
+                            [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Ev2},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
+                            for n=1:length(PSR2)
+                                Ev2Spikes(n,1)=sum(PSR2{1,n}>(binstart));
+                            end
+                            
+                            %pick which trials get used for decoding
+                            SetupTrials2=cat(1,ones(Ev2s,1),zeros(length(PSR2)-Ev2s,1));
+                            Trials2=SetupTrials2==1;%(SetupTrials2(randperm(length(SetupTrials2)))==1); %randomize trials (or not)
+                            
+                            %put spikes from those trials in a matrix
+                            AllSpikes((Ev1s+1):(Ev1s+Ev2s),NN)=Ev2Spikes(Trials2,1);
+                            NN=NN+1; %neuron counter
+                            AllNN=AllNN+1;
+                        else AllNN=AllNN+1;
+                        end
+                    end
+                end
+                
+                
+                
+                %Setup spikes matrix for decoding
+                DecodeSpikes=AllSpikes(:,Sel);
+                
+                %find neurons with too few spikes
+                
+                for t=1:NumNeurons(v)
+                    if sum(DecodeSpikes(1:Ev1s,t))<3 || sum(DecodeSpikes((1+Ev1s):(Ev1s+Ev2s),t))<3
+                        LowHz(1,t)=1;
+                    end
+                end
+                
+                %remove neurons with too few spikes
+%                 if sum(LowHz)>0
+%                     Sel2=LowHz<1;
+%                     DecodeSpikes=DecodeSpikes(:,Sel2);
+%                 end
+                
+                %set up validation result matrices
+                CVacc = NaN(folds,1);
+                CVaccSh = NaN(folds,1);
+                
+                %normal model
+                Partitions = cvpartition(DecodeRs,'KFold',folds);
+                for r = 1:folds
+                    
+                    SVMModel = fitcdiscr(DecodeSpikes(Partitions.training(r),:),DecodeRs(Partitions.training(r)),'discrimType','diaglinear');
+                    prediction = predict(SVMModel,DecodeSpikes(Partitions.test(r),:));
+                    actual = DecodeRs(Partitions.test(r));
+                    correct = prediction - actual;
+                    CVacc(r) = sum(correct==0) / length(correct);
+                    
+                end
+                
+                %shuffled model
+                for q=1:shuffs
+                    DecodeRsSh=DecodeRs(randperm(length(DecodeRs)));
+                    PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
+                    for s = 1:folds
+                        SVMModelSh = fitcdiscr(DecodeSpikes(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)),'discrimType','diaglinear');
+                        predictionSh = predict(SVMModelSh,DecodeSpikes(PartitionsSh.test(s),:));
+                        actualSh = DecodeRsSh(PartitionsSh.test(s));
+                        correctSh = predictionSh - actualSh;
+                        CVaccSh(s) = sum(correctSh==0) / length(correctSh);
+                    end
+                    AccShuff(q,1) = nanmean(CVaccSh);
+                end
+                
+                PoolMdlAcc{v,1}(u,l) = nanmean(CVacc);
+                PoolMdlAccShuff{v,1}(u,l) = nanmean(AccShuff);
+                
+                fprintf(['Condition #' num2str(v) ', Rep #' num2str(u) ', Bin #' num2str(l) '\n']);
+            end %bins
+        end
+    end
+    
+    
+    
+    if PavInst==1
+        PoolMdlAccDSSucrose=PoolMdlAcc;
+        PoolMdlAccDSSucroseShuff=PoolMdlAccShuff;
+        save('PoolMdlAccDSSucrose.mat','PoolMdlAccDSSucrose')
+        save('PoolMdlAccDSSucroseShuff.mat','PoolMdlAccDSSucroseShuff')
+    else if PavInst==2
+            PoolMdlAccPavSucrose=PoolMdlAcc;
+            PoolMdlAccPavSucroseShuff=PoolMdlAccShuff;
+            save('PoolMdlAccPavSucrose.mat','PoolMdlAccPavSucrose')
+            save('PoolMdlAccPavSucroseShuff.mat','PoolMdlAccPavSucroseShuff')
+        else if PavInst==3
+                PoolMdlAccDSAlcohol=PoolMdlAcc;
+                PoolMdlAccDSAlcoholShuff=PoolMdlAccShuff;
+                save('PoolMdlAccDSAlcohol.mat','PoolMdlAccDSAlcohol')
+                save('PoolMdlAccDSAlcoholShuff.mat','PoolMdlAccDSAlcoholShuff')
+            else if PavInst==4
+                    PoolMdlAccPavAlcohol=PoolMdlAcc;
+                    PoolMdlAccPavAlcoholShuff=PoolMdlAccShuff;
+                    save('PoolMdlAccPavAlcohol.mat','PoolMdlAccPavAlcohol')
+                    save('PoolMdlAccPavAlcoholShuff.mat','PoolMdlAccPavAlcoholShuff')
+                    
+                end
+            end
+        end
+    end
+    
+end
+
+%%
+% 
+% for PavInst=1:4;
+%     if PavInst==1
+%         PoolMdlAccDSSucroseQuadratic=PoolMdlAccDSSucrose;
+%         PoolMdlAccDSSucroseShuffQuadratic=PoolMdlAccDSSucroseShuff;
+%         save('PoolMdlAccDSSucroseQuadratic.mat','PoolMdlAccDSSucroseQuadratic')
+%         save('PoolMdlAccDSSucroseShuffQuadratic.mat','PoolMdlAccDSSucroseShuffQuadratic')
+%     else if PavInst==2
+%             PoolMdlAccPavSucroseQuadratic=PoolMdlAccPavSucrose;
+%             PoolMdlAccPavSucroseShuffQuadratic=PoolMdlAccPavSucroseShuff;
+%             save('PoolMdlAccPavSucroseQuadratic.mat','PoolMdlAccPavSucroseQuadratic')
+%             save('PoolMdlAccPavSucroseShuffQuadratic.mat','PoolMdlAccPavSucroseShuffQuadratic')
+%         else if PavInst==3
+%                 PoolMdlAccDSAlcoholQuadratic=PoolMdlAccDSAlcohol;
+%                 PoolMdlAccDSAlcoholShuffQuadratic=PoolMdlAccDSAlcoholShuff;
+%                 save('PoolMdlAccDSAlcoholQuadratic.mat','PoolMdlAccDSAlcoholQuadratic')
+%                 save('PoolMdlAccDSAlcoholShuffQuadratic.mat','PoolMdlAccDSAlcoholQuadratic')
+%             else if PavInst==4
+%                     PoolMdlAccPavAlcoholQuadratic=PoolMdlAccPavAlcohol;
+%                     PoolMdlAccPavAlcoholShuffQuadratic=PoolMdlAccPavAlcoholShuff;
+%                     save('PoolMdlAccPavAlcoholQuadratic.mat','PoolMdlAccPavAlcoholQuadratic')
+%                     save('PoolMdlAccPavAlcoholShuffQuadratic.mat','PoolMdlAccPavAlcoholShuffQuadratic')
+%                     
+%                 end
+%             end
+%         end
+%     end
+% end

+ 453 - 0
AnalysisAndFigures/o_Figure7_DecodingAlcoholExposed.m

@@ -0,0 +1,453 @@
+%Plotting decoding analysis
+clear all;
+load('UnitMdlDSSucrose.mat')
+load('UnitMdlDSShuffSucrose.mat')
+load('UnitMdlPavSucrose.mat')
+load('UnitMdlPavShuffSucrose.mat')
+if ~exist('RPPSall'),load('RPPSall.mat'); end
+if ~exist('RPSall'),load('RPSall.mat'); end
+load('UnitMdlDSAlcohol.mat')
+load('UnitMdlDSShuffAlcohol.mat')
+load('UnitMdlPavAlcohol.mat')
+load('UnitMdlPavShuffAlcohol.mat')
+if ~exist('RvppasALL'),load('RvppasALL.mat'); end
+if ~exist('RvpppasALL'),load('RvpppasALL.mat'); end
+
+sucrose=[1  0.6  0.1];
+maltodextrin=[.9  0.3  .9];
+water=[0.00  0.75  0.75];
+total=[0.3  0.1  0.8];
+DS=[0    0.4470    0.7410];
+Pav=[0.8500    0.3250    0.0980];
+DSAlcohol=[51/255 160/255 44/255];%[0.4940    0.1840    0.5560];
+PavAlcohol=[180/255 41/255 230/255];%[0.9290    0.6940    0.1250];
+
+DSShuff=[0 0 0];
+PavShuff=[0 0 0];
+
+xaxis=linspace(-0.75,3.15,13);
+
+selection1=RPPSall.Structure==10 & RPPSall.CSPlusRatio>=.5 & (RPPSall.CSPlusRatio./(RPPSall.CSPlusRatio+RPPSall.CSMinusRatio))>=.7;% & RPPSall.Ev(1).RespDir==1;
+selection2=RPSall.Structure==10 & RPSall.CSPlusRatio>=.5 & (RPSall.CSPlusRatio./(RPSall.CSPlusRatio+RPSall.CSMinusRatio))>=.7;% & RPSall.Ev(1).RespDir==1;
+selection3=RvpppasALL.Structure==10 & RvpppasALL.CSPlusRatio>=.5 & (RvpppasALL.CSPlusRatio./(RvpppasALL.CSPlusRatio+RvpppasALL.CSMinusRatio))>=.7;% & RvpppasALL.Ev(1).RespDir==1;
+selection4=RvppasALL.Structure==10 & RvppasALL.CSPlusRatio>=.5 & (RvppasALL.CSPlusRatio./(RvppasALL.CSPlusRatio+RvppasALL.CSMinusRatio))>=.7;% & RvppasALL.Ev(1).RespDir==1;
+
+%Figure code
+
+figure;
+
+%get average accuracy for each bin
+for i = 1:13
+    AvgAccDSSucrose(1,i)=nanmean(UnitMdlAccDSSucrose(selection1,i)); %average accuracy DS
+    SEMAccDSSucrose(1,i)=nanste(UnitMdlAccDSSucrose(selection1,i),1); %SEM
+    AvgAccDSShuffSucrose(1,i)=nanmean(UnitMdlAccDSShuffSucrose(selection1,i)); %average accuracy DSShuff
+    SEMAccDSShuffSucrose(1,i)=nanste(UnitMdlAccDSShuffSucrose(selection1,i),1); %SEM
+    AvgAccPavSucrose(1,i)=nanmean(UnitMdlAccPavSucrose(selection2,i)); %average accuracy DSShuff
+    SEMAccPavSucrose(1,i)=nanste(UnitMdlAccPavSucrose(selection2,i),1); %SEM
+    AvgAccPavShuffSucrose(1,i)=nanmean(UnitMdlAccPavShuffSucrose(selection2,i)); %average accuracy DSShuff
+    SEMAccPavShuffSucrose(1,i)=nanste(UnitMdlAccPavShuffSucrose(selection2,i),1); %SEM
+    
+    DSSigSucrose(1,i) = ttest2(UnitMdlAccDSSucrose(:,i),UnitMdlAccDSShuffSucrose(selection1,i),'Vartype','unequal'); %welch's t-test
+    PavSigSucrose(1,i) = ttest2(UnitMdlAccPavSucrose(:,i),UnitMdlAccPavShuffSucrose(selection2,i),'Vartype','unequal'); %welch's t-test
+end
+
+%get normalized accuracy increase over shuffled for each neuron, then take
+%average and find SEM
+for i = 1:13
+    for j = 1:length(UnitMdlAccDSSucrose)
+        NormAccDSSucrose(j,i) = UnitMdlAccDSSucrose(j,i)-nanmean(nanmean(UnitMdlAccDSShuffSucrose(:,i)));
+    end
+    for k = 1:length(UnitMdlAccPavSucrose)
+        NormAccPavSucrose(k,i) = UnitMdlAccPavSucrose(k,i)-nanmean(nanmean(UnitMdlAccPavShuffSucrose(:,i)));
+    end
+    
+    AvgNormDSSucrose(1,i)=nanmean(NormAccDSSucrose(selection1,i)); %average norm accuracy DS
+    SEMNormDSSucrose(1,i)=nanste(NormAccDSSucrose(selection1,i),1); %SEM
+    AvgNormPavSucrose(1,i)=nanmean(NormAccPavSucrose(selection2,i)); %average norm accuracy DS
+    SEMNormPavSucrose(1,i)=nanste(NormAccPavSucrose(selection2,i),1); %SEM
+    
+    %NormSig(1,i) = ttest2(NormAccDS(:,i),NormAccPav(:,i),'Vartype','unequal'); %welch's t-test
+    [~,NormSig(1,i)] = ranksum(NormAccDSSucrose(selection1,i),NormAccPavSucrose(selection2,i)); %wilcoxon's ranksum test
+end
+
+%prepare shading
+upSEMDS1=AvgAccDSSucrose+SEMAccDSSucrose;
+downSEMDS1=AvgAccDSSucrose-SEMAccDSSucrose;
+upSEMPav1=AvgAccPavSucrose+SEMAccPavSucrose;
+downSEMPav1=AvgAccPavSucrose-SEMAccPavSucrose;
+upSEMDSShuff1=AvgAccDSShuffSucrose+SEMAccDSShuffSucrose;
+downSEMDSShuff1=AvgAccDSShuffSucrose-SEMAccDSShuffSucrose;
+upSEMPavShuff1=AvgAccPavShuffSucrose+SEMAccPavShuffSucrose;
+downSEMPavShuff1=AvgAccPavShuffSucrose-SEMAccPavShuffSucrose;
+upSEMNormDS1=AvgNormDSSucrose+SEMNormDSSucrose;
+downSEMNormDS1=AvgNormDSSucrose-SEMNormDSSucrose;
+upSEMNormPav1=AvgNormPavSucrose+SEMNormPavSucrose;
+downSEMNormPav1=AvgNormPavSucrose-SEMNormPavSucrose;
+
+% %plotting decoder accuracy over time
+% subplot(3,4,1); %accumbens
+% hold on;
+% plot(xaxis,AvgAccDSSucrose(1:13),'Color', DS,'linewidth',3);
+% plot(xaxis,AvgAccDSShuffSucrose(1:13),'Color', DSShuff,'linewidth',3);
+% %plot(xaxis,DSSigSucrose-0.51,'*','Color', 'k');
+% patch([xaxis,xaxis(end:-1:1)],[upSEMDS,downSEMDS(end:-1:1)],DS,'EdgeColor','none');alpha(0.5);
+% patch([xaxis,xaxis(end:-1:1)],[upSEMDSShuff,downSEMDSShuff(end:-1:1)],DSShuff,'EdgeColor','none');alpha(0.5);
+% xlabel('Seconds post cue');
+% legend('Sucrose DS','Shuffled','Location','best');
+% %title('Sucrose DS single unit decoding accuracy');
+% ylabel('Single Unit Decoding Accuracy');
+% axis([-.75 3.15 0.45 0.7]);
+% 
+% subplot(3,4,2); %Pav
+% hold on;
+% plot(xaxis,AvgAccPavSucrose(1:13),'Color', Pav,'linewidth',3);
+% plot(xaxis,AvgAccPavShuffSucrose(1:13),'Color', PavShuff,'linewidth',3);
+% %plot(xaxis,PavSig-0.51,'*','Color', 'k');
+% patch([xaxis,xaxis(end:-1:1)],[upSEMPav,downSEMPav(end:-1:1)],Pav,'EdgeColor','none');alpha(0.5);
+% patch([xaxis,xaxis(end:-1:1)],[upSEMPavShuff,downSEMPavShuff(end:-1:1)],PavShuff,'EdgeColor','none');alpha(0.5);
+% xlabel('Seconds post cue');
+% %title('Sucrose CS+ single unit decoding accuracy');
+% legend('Sucrose CS+','Shuffled','Location','northeast');
+% axis([-.75 3.15 0.45 0.7]);
+
+%get average accuracy for each bin
+for i = 1:13
+    AvgAccDSAlcohol(1,i)=nanmean(UnitMdlAccDSAlcohol(selection3,i)); %average accuracy DS
+    SEMAccDSAlcohol(1,i)=nanste(UnitMdlAccDSAlcohol(selection3,i),1); %SEM
+    AvgAccDSShuffAlcohol(1,i)=nanmean(UnitMdlAccDSShuffAlcohol(selection3,i)); %average accuracy DSShuff
+    SEMAccDSShuffAlcohol(1,i)=nanste(UnitMdlAccDSShuffAlcohol(selection3,i),1); %SEM
+    AvgAccPavAlcohol(1,i)=nanmean(UnitMdlAccPavAlcohol(selection4,i)); %average accuracy DSShuff
+    SEMAccPavAlcohol(1,i)=nanste(UnitMdlAccPavAlcohol(selection4,i),1); %SEM
+    AvgAccPavShuffAlcohol(1,i)=nanmean(UnitMdlAccPavShuffAlcohol(selection4,i)); %average accuracy DSShuff
+    SEMAccPavShuffAlcohol(1,i)=nanste(UnitMdlAccPavShuffAlcohol(selection4,i),1); %SEM
+    
+    DSSigAlcohol(1,i) = ttest2(UnitMdlAccDSAlcohol(:,i),UnitMdlAccDSShuffAlcohol(selection3,i),'Vartype','unequal'); %welch's t-test
+    PavSigAlcohol(1,i) = ttest2(UnitMdlAccPavAlcohol(:,i),UnitMdlAccPavShuffAlcohol(selection4,i),'Vartype','unequal'); %welch's t-test
+end
+
+%get normalized accuracy increase over shuffled for each neuron, then take
+%average and find SEM
+for i = 1:13
+    for j = 1:length(UnitMdlAccDSAlcohol)
+        NormAccDSAlcohol(j,i) = UnitMdlAccDSAlcohol(j,i)-nanmean(nanmean(UnitMdlAccDSShuffAlcohol(:,i)));
+    end
+    for k = 1:length(UnitMdlAccPavAlcohol)
+        NormAccPavAlcohol(k,i) = UnitMdlAccPavAlcohol(k,i)-nanmean(nanmean(UnitMdlAccPavShuffAlcohol(:,i)));
+    end
+    
+    AvgNormDSAlcohol(1,i)=nanmean(NormAccDSAlcohol(selection3,i)); %average norm accuracy DS
+    SEMNormDSAlcohol(1,i)=nanste(NormAccDSAlcohol(selection3,i),1); %SEM
+    AvgNormPavAlcohol(1,i)=nanmean(NormAccPavAlcohol(selection4,i)); %average norm accuracy DS
+    SEMNormPavAlcohol(1,i)=nanste(NormAccPavAlcohol(selection4,i),1); %SEM
+    
+    %NormSig(1,i) = ttest2(NormAccDS(:,i),NormAccPav(:,i),'Vartype','unequal'); %welch's t-test
+    [~,NormSig(1,i)] = ranksum(NormAccDSAlcohol(selection3,i),NormAccPavAlcohol(selection4,i)); %wilcoxon's ranksum test
+end
+
+%prepare shading
+upSEMDS=AvgAccDSAlcohol+SEMAccDSAlcohol;
+downSEMDS=AvgAccDSAlcohol-SEMAccDSAlcohol;
+upSEMPav=AvgAccPavAlcohol+SEMAccPavAlcohol;
+downSEMPav=AvgAccPavAlcohol-SEMAccPavAlcohol;
+upSEMDSShuff=AvgAccDSShuffAlcohol+SEMAccDSShuffAlcohol;
+downSEMDSShuff=AvgAccDSShuffAlcohol-SEMAccDSShuffAlcohol;
+upSEMPavShuff=AvgAccPavShuffAlcohol+SEMAccPavShuffAlcohol;
+downSEMPavShuff=AvgAccPavShuffAlcohol-SEMAccPavShuffAlcohol;
+upSEMNormDS=AvgNormDSAlcohol+SEMNormDSAlcohol;
+downSEMNormDS=AvgNormDSAlcohol-SEMNormDSAlcohol;
+upSEMNormPav=AvgNormPavAlcohol+SEMNormPavAlcohol;
+downSEMNormPav=AvgNormPavAlcohol-SEMNormPavAlcohol;
+
+% %plotting decoder accuracy over time
+% subplot(3,4,5); %accumbens
+% hold on;
+% plot(xaxis,AvgAccDSAlcohol(1:13),'Color', DSAlcohol,'linewidth',3);
+% plot(xaxis,AvgAccDSShuffAlcohol(1:13),'Color', DSShuff,'linewidth',3);
+% %plot(xaxis,DSSigAlcohol-0.51,'*','Color', 'k');
+% patch([xaxis,xaxis(end:-1:1)],[upSEMDS,downSEMDS(end:-1:1)],DSAlcohol,'EdgeColor','none');alpha(0.5);
+% patch([xaxis,xaxis(end:-1:1)],[upSEMDSShuff,downSEMDSShuff(end:-1:1)],DSShuff,'EdgeColor','none');alpha(0.5);
+% xlabel('Seconds post cue');
+% legend('Suc(alc) DS','Shuffled','Location','northeast');
+% %title('Alcohol DS single unit decoding accuracy');
+% axis([-.75 3.15 0.45 0.7]);
+% ylabel('Single Unit Decoding Accuracy');
+% 
+% subplot(3,4,6); %Pav
+% hold on;
+% plot(xaxis,AvgAccPavAlcohol(1:13),'Color', PavAlcohol,'linewidth',3);
+% plot(xaxis,AvgAccPavShuffAlcohol(1:13),'Color', PavShuff,'linewidth',3);
+% %plot(xaxis,PavSig-0.51,'*','Color', 'k');
+% patch([xaxis,xaxis(end:-1:1)],[upSEMPav,downSEMPav(end:-1:1)],PavAlcohol,'EdgeColor','none');alpha(0.5);
+% patch([xaxis,xaxis(end:-1:1)],[upSEMPavShuff,downSEMPavShuff(end:-1:1)],PavShuff,'EdgeColor','none');alpha(0.5);
+% xlabel('Seconds post cue');
+% %title('Alcohol CS+ single unit decoding accuracy');
+% legend('Suc(alc) CS+','Shuffled','Location','northeast');
+% axis([-.75 3.15 0.45 0.7]);
+
+% combined single unit decoding over time
+subplot(2,2,1);
+hold on;
+plot(xaxis,AvgAccDSAlcohol(1:13),'Marker','o','Color', DSAlcohol,'linewidth',1);
+plot(xaxis,AvgAccPavAlcohol(1:13),'Marker','o','Color', PavAlcohol,'linewidth',1);
+plot(xaxis,AvgAccDSShuffAlcohol(1:13),'Marker','o','Color', DSShuff,'linewidth',1);
+plot(xaxis,AvgAccPavShuffAlcohol(1:13),'Marker','o','Color', PavShuff,'linewidth',1);
+%plot(xaxis,DSSigAlcohol-0.51,'*','Color', 'k');
+patch([0 0.3 0.3 0],[0 0 1 1],[0.5 0.5 0.5],'EdgeColor','none');alpha(0.5);
+patch([xaxis,xaxis(end:-1:1)],[upSEMDS,downSEMDS(end:-1:1)],DSAlcohol,'EdgeColor','none');alpha(0.5);
+patch([xaxis,xaxis(end:-1:1)],[upSEMDSShuff,downSEMDSShuff(end:-1:1)],DSShuff,'EdgeColor','none');alpha(0.5);
+patch([xaxis,xaxis(end:-1:1)],[upSEMPav,downSEMPav(end:-1:1)],PavAlcohol,'EdgeColor','none');alpha(0.5);
+patch([xaxis,xaxis(end:-1:1)],[upSEMPavShuff,downSEMPavShuff(end:-1:1)],PavShuff,'EdgeColor','none');alpha(0.5);
+
+
+%this adds all 4 plots
+% plot(xaxis,AvgAccDSSucrose(1:13),'Color', DS,'linewidth',1);
+% plot(xaxis,AvgAccDSShuffSucrose(1:13),'Color', DSShuff,'linewidth',1);
+% patch([xaxis,xaxis(end:-1:1)],[upSEMDS1,downSEMDS1(end:-1:1)],DS,'EdgeColor','none');alpha(0.5);
+% patch([xaxis,xaxis(end:-1:1)],[upSEMDSShuff1,downSEMDSShuff1(end:-1:1)],DSShuff,'EdgeColor','none');alpha(0.5);
+% plot(xaxis,AvgAccPavSucrose(1:13),'Color', Pav,'linewidth',1);
+% plot(xaxis,AvgAccPavShuffSucrose(1:13),'Color', PavShuff,'linewidth',1);
+% patch([xaxis,xaxis(end:-1:1)],[upSEMPav1,downSEMPav1(end:-1:1)],Pav,'EdgeColor','none');alpha(0.5);
+% patch([xaxis,xaxis(end:-1:1)],[upSEMPavShuff1,downSEMPavShuff1(end:-1:1)],PavShuff,'EdgeColor','none');alpha(0.5);
+
+xlabel('Seconds post cue');
+legend('Suc(alc) DS','Suc(alc) CS+','Shuffled','Location','northeast');
+%title('Alcohol DS single unit decoding accuracy');
+axis([-.75 3.15 0.45 0.7]);
+ylabel('Single Unit Decoding Accuracy');
+
+%plotting cdf plots for DS at time of max predictive power
+subplot(2,2,2);
+DSbin=4;
+[cdfDS,xDS] = ecdf(UnitMdlAccDSAlcohol(selection3,DSbin));
+[cdfPav,xPav] = ecdf(UnitMdlAccPavAlcohol(selection4,DSbin));
+[cdfDSsuc,xDSsuc] = ecdf(UnitMdlAccDSSucrose(selection1,DSbin));
+[cdfPavsuc,xPavsuc] = ecdf(UnitMdlAccPavSucrose(selection2,DSbin));
+%subplot(4,3,[3 6]);
+hold on;
+plot(xDSsuc,cdfDSsuc,'Color',DS,'linewidth',1);
+plot(xPavsuc,cdfPavsuc,'Color',Pav,'linewidth',1);
+plot(xDS,cdfDS,'Color',DSAlcohol,'linewidth',1);
+plot(xPav,cdfPav,'Color',PavAlcohol,'linewidth',1);
+plot([.5 .5], [0 1],'LineStyle',':','Color','k');
+plot([AvgAccDSSucrose(1,4) AvgAccDSSucrose(1,4)], [0 1],'Color',DS,'linewidth',1);
+plot([AvgAccDSAlcohol(1,4) AvgAccDSAlcohol(1,4)], [0 1],'Color',DSAlcohol,'linewidth',1);
+plot([AvgAccPavSucrose(1,4) AvgAccPavSucrose(1,4)], [0 1],'Color',Pav,'linewidth',1);
+plot([AvgAccPavAlcohol(1,4) AvgAccPavAlcohol(1,4)], [0 1],'Color',PavAlcohol,'linewidth',1);
+axis([0.25 1 0 1]);
+title(['Cumulative Distribution 0 to .3s']);
+xlabel('Decoding Accuracy')
+ylabel('Proportion of Population')
+legend('Sucrose DS','Sucrose CS+','Suc(alc) DS','Suc(alc) CS+','Location','best');
+
+
+%Plotting ensemble decoding
+
+clear all;
+DS=[0    0.4470    0.7410];
+Pav=[0.8500    0.3250    0.0980];
+DSAlcohol=[51/255 160/255 44/255];%[0.4940    0.1840    0.5560];
+PavAlcohol=[180/255 41/255 230/255];%[0.9290    0.6940    0.1250];
+NumNeurons=[1 5 10 25 50 100];
+
+load('PoolMdlAccDSSucrose.mat','PoolMdlAccDSSucrose')
+load('PoolMdlAccDSSucroseShuff.mat','PoolMdlAccDSSucroseShuff')
+load('PoolMdlAccPavSucrose.mat','PoolMdlAccPavSucrose')
+load('PoolMdlAccPavSucroseShuff.mat','PoolMdlAccPavSucroseShuff')
+load('PoolMdlAccDSAlcohol.mat','PoolMdlAccDSAlcohol')
+load('PoolMdlAccDSAlcoholShuff.mat','PoolMdlAccDSAlcoholShuff')
+load('PoolMdlAccPavAlcohol.mat','PoolMdlAccPavAlcohol')
+load('PoolMdlAccPavAlcoholShuff.mat','PoolMdlAccPavAlcoholShuff')
+
+
+for i=1:6    
+AvgAccDSSucrose(i)=nanmean(PoolMdlAccDSSucrose{i,1}); %average accuracy
+SEMAccDSSucrose(i)=nanste(PoolMdlAccDSSucrose{i,1},1); %SEM
+
+AvgAccPavSucrose(i)=nanmean(PoolMdlAccPavSucrose{i,1}); %average accuracy
+SEMAccPavSucrose(i)=nanste(PoolMdlAccPavSucrose{i,1},1); %SEM
+
+AvgAccDSAlcohol(i)=nanmean(PoolMdlAccDSAlcohol{i,1}); %average accuracy
+SEMAccDSAlcohol(i)=nanste(PoolMdlAccDSAlcohol{i,1},1); %SEM
+
+AvgAccPavAlcohol(i)=nanmean(PoolMdlAccPavAlcohol{i,1}); %average accuracy
+SEMAccPavAlcohol(i)=nanste(PoolMdlAccPavAlcohol{i,1},1); %SEM
+
+AvgAccDSSucroseShuff(i)=nanmean(PoolMdlAccDSSucroseShuff{i,1}); %average accuracy
+SEMAccDSSucroseShuff(i)=nanste(PoolMdlAccDSSucroseShuff{i,1},1); %SEM
+
+AvgAccPavSucroseShuff(i)=nanmean(PoolMdlAccPavSucroseShuff{i,1}); %average accuracy
+SEMAccPavSucroseShuff(i)=nanste(PoolMdlAccPavSucroseShuff{i,1},1); %SEM
+
+AvgAccDSAlcoholShuff(i)=nanmean(PoolMdlAccDSAlcoholShuff{i,1}); %average accuracy
+SEMAccDSAlcoholShuff(i)=nanste(PoolMdlAccDSAlcoholShuff{i,1},1); %SEM
+
+AvgAccPavAlcoholShuff(i)=nanmean(PoolMdlAccPavAlcoholShuff{i,1}); %average accuracy
+SEMAccPavAlcoholShuff(i)=nanste(PoolMdlAccPavAlcoholShuff{i,1},1); %SEM
+end
+
+
+
+subplot(2,2,3); %ensemble decoding results
+errorbar(NumNeurons, AvgAccDSSucrose,SEMAccDSSucrose,'Color',DS,'linewidth',1);
+hold on;
+errorbar(NumNeurons, AvgAccPavSucrose,SEMAccPavSucrose, 'Color',Pav,'linewidth',1);
+errorbar(NumNeurons, AvgAccDSAlcohol,SEMAccDSAlcohol, 'Color',DSAlcohol,'linewidth',1);
+errorbar(NumNeurons, AvgAccPavAlcohol,SEMAccPavAlcohol, 'Color',PavAlcohol,'linewidth',1);
+% errorbar(NumNeurons, AvgAccDSSucroseShuff,SEMAccDSSucroseShuff, 'Color','k');
+% errorbar(NumNeurons, AvgAccPavSucroseShuff,SEMAccPavSucroseShuff, 'Color','k');
+% errorbar(NumNeurons, AvgAccDSAlcoholShuff,SEMAccDSAlcoholShuff, 'Color','k');
+% errorbar(NumNeurons, AvgAccPavAlcoholShuff,SEMAccPavAlcoholShuff, 'Color','k');
+legend('Sucrose DS','Sucrose CS+','Suc(alc) DS','Suc(alc) CS+','Location','best');
+plot([-5 105], [.5 .5],'LineStyle',':','Color','k');
+xlabel('Neurons used');
+ylabel('Mean Accuracy (+/- SEM)');
+xticks([1 5 10 25 50 100])
+%xticklabels({'-3\pi','-2\pi','-\pi','0','\pi','2\pi','3\pi'})
+axis([-5 105 0.4 1]);
+
+%make matrices with bins at all levels
+A=[];
+B=[];
+C=[];
+D=[];
+A1=[];
+B1=[];
+C1=[];
+D1=[];
+
+for v=1:length(PoolMdlAccDSSucrose)
+A=cat(1,A,PoolMdlAccDSSucrose{v,1});
+B=cat(1,B,PoolMdlAccPavSucrose{v,1});
+C=cat(1,C,PoolMdlAccDSAlcohol{v,1});
+D=cat(1,D,PoolMdlAccPavAlcohol{v,1});
+A1=cat(1,A1,PoolMdlAccDSSucroseShuff{v,1});
+B1=cat(1,B1,PoolMdlAccPavSucroseShuff{v,1});
+C1=cat(1,C1,PoolMdlAccDSAlcoholShuff{v,1});
+D1=cat(1,D1,PoolMdlAccPavAlcoholShuff{v,1});
+end
+
+subplot(2,2,4); %summary of accuracy
+errorbar(1,nanmean(A),nanste(A,1),'o','Color',DS,'linewidth',1);
+hold on;
+errorbar(2,nanmean(C),nanste(C,1),'o', 'Color',DSAlcohol,'linewidth',1);
+errorbar(3,nanmean(B),nanste(B,1),'o', 'Color',Pav,'linewidth',1);
+errorbar(4,nanmean(D),nanste(D,1),'o', 'Color',PavAlcohol,'linewidth',1);
+plot([-5 105], [.5 .5],'LineStyle',':','Color','k');
+plot([1 2], [0.93 0.93],'Color','k','linewidth',1);
+plot([3 4], [0.86 0.86],'Color','k','linewidth',1);
+plot(1.5,.95,'*','Color','k');
+plot(3.5,.88,'*','Color','k');
+axis([0.5 4.5 0.4 1]);
+set(gca,'xtick',[])
+ylabel('Mean Accuracy (+/- SEM)');
+%xlabel('DS/Suc DS/Suc(a) CS+/Suc CS+/Suc(a)');
+legend('Sucrose DS','Suc(alc) DS','Sucrose CS+','Suc(alc) CS+','Location','best');
+
+%% STATS
+
+%Plotting single All decoding
+bins = 1; %number of bins
+binint = 0.3; %spacing of bins
+binstart = 0; %start time of first bin relative to event
+Dura=[0 0.3]; %size of bin used for decoding:
+NumNeurons=[1 5 10 25 50 100]; %matrix of how many neurons to use on each iteration
+repetitions=20; %how many times to run the analysis
+
+%ANOVA for effect of ensemble size, task and reward on accuracy 0-.3 sec
+size=[];
+for i=1:length(NumNeurons)
+    size=cat(1,size,i*ones(repetitions,1));
+end
+size2=cat(1,size,size,size,size);
+%shuffd=cat(1,zeros(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1));
+reward=cat(1,ones(length(A)+length(B),1),zeros(length(C)+length(D),1));
+task=cat(1,ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1));
+[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D),{reward,task,size2},'varnames',{'reward','task','ens size'},'model','full');
+%multcompare(AccStats,'Dimension',[1:3])
+
+rewardtask=cat(1,zeros(length(A),1),ones(length(B),1),2*ones(length(C),1),3*ones(length(D),1));
+[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D),{rewardtask,size2},'varnames',{'reward and task','ens size'},'model','full');
+
+
+size=[];
+for i=1:length(NumNeurons)
+    size=cat(1,size,i*ones(repetitions,1));
+end
+size1=cat(1,size,size,size,size,size,size,size,size);
+shuffd=cat(1,zeros(length(A)+length(B)+length(C)+length(D),1),ones(length(A)+length(B)+length(C)+length(D),1));
+reward=cat(1,ones(length(A)+length(B),1),zeros(length(C)+length(D),1),ones(length(A)+length(B),1),zeros(length(C)+length(D),1));
+task=cat(1,ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1),ones(length(A),1),zeros(length(B),1),ones(length(C),1),zeros(length(D),1));
+[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D,A1,B1,C1,D1),{reward,task,size1,shuffd},'varnames',{'reward','task','ens size','shuffled vs real'},'model','full');
+
+
+rewardtask=cat(1,zeros(length(A),1),ones(length(B),1),2*ones(length(C),1),3*ones(length(D),1));
+[~,AccAnova,AccStats]=anovan(cat(1,A,B,C,D),{rewardtask,size2},'varnames',{'reward and task','ens size'},'model','full');
+
+
+%without shuff vs true
+% size3=cat(1,size,size);
+% region3=cat(1,ones(length(A),1),zeros(length(C),1));
+% [~,AccAnovaNoShuff,~]=anovan(cat(1,A,C),{region3,size3},'varnames',{'region','ens size'},'display','off','model','full');
+
+%% Sinlge Unit Stats
+
+if ~exist('RPPSall'),load('RPPSall.mat'); end
+if ~exist('RPSall'),load('RPSall.mat'); end
+load('UnitMdlDSAlcohol.mat')
+load('UnitMdlDSShuffAlcohol.mat')
+load('UnitMdlPavAlcohol.mat')
+load('UnitMdlPavShuffAlcohol.mat')
+if ~exist('RvppasALL'),load('RvppasALL.mat'); end
+if ~exist('RvpppasALL'),load('RvpppasALL.mat'); end
+
+sucrose=[1  0.6  0.1];
+maltodextrin=[.9  0.3  .9];
+water=[0.00  0.75  0.75];
+total=[0.3  0.1  0.8];
+DS=[0    0.4470    0.7410];
+Pav=[0.8500    0.3250    0.0980];
+DSAlcohol=[0.4940    0.1840    0.5560];
+PavAlcohol=[0.9290    0.6940    0.1250];
+DSShuff=[0 0 0];
+PavShuff=[0 0 0];
+
+xaxis=linspace(-0.75,3.15,13);
+
+selection1=RPPSall.Structure==10 & RPPSall.CSPlusRatio>=.5 & (RPPSall.CSPlusRatio./(RPPSall.CSPlusRatio+RPPSall.CSMinusRatio))>=.7;
+selection2=RPSall.Structure==10 & RPSall.CSPlusRatio>=.5 & (RPSall.CSPlusRatio./(RPSall.CSPlusRatio+RPSall.CSMinusRatio))>=.7;
+selection3=RvpppasALL.Structure==10 & RvpppasALL.CSPlusRatio>=.5 & (RvpppasALL.CSPlusRatio./(RvpppasALL.CSPlusRatio+RvpppasALL.CSMinusRatio))>=.7;
+selection4=RvppasALL.Structure==10 & RvppasALL.CSPlusRatio>=.5 & (RvppasALL.CSPlusRatio./(RvppasALL.CSPlusRatio+RvppasALL.CSMinusRatio))>=.7;
+
+%LME for the 2 alc exposed
+task=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccDSAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1));
+shuffd=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),zeros(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccPavAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1));
+rat=cat(1,RvpppasALL.Rat,RvpppasALL.Rat,RvppasALL.Rat,RvppasALL.Rat);
+tbl=table(cat(1,UnitMdlAccDSAlcohol(:,4),UnitMdlAccDSShuffAlcohol(:,4),UnitMdlAccPavAlcohol(:,4),UnitMdlAccPavShuffAlcohol(:,4)),task,shuffd,rat,'variablenames',{'Accuracy','Task','Shuffled','Subject'});
+lme=fitlme(tbl,'Accuracy~Task+Shuffled+(Task:Shuffled)+(1|Subject)');%+(1|Subject:Task:Shuffled)');
+
+%all 4
+load('UnitMdlDSSucrose.mat')
+load('UnitMdlDSShuffSucrose.mat')
+load('UnitMdlPavSucrose.mat')
+load('UnitMdlPavShuffSucrose.mat')
+if ~exist('RPPSall'),load('RPPSall.mat'); end
+if ~exist('RPSall'),load('RPSall.mat'); end
+task=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccDSAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1),ones(length(UnitMdlAccDSSucrose(:,4)),1),ones(length(UnitMdlAccDSSucrose(:,4)),1),zeros(length(UnitMdlAccPavSucrose(:,4)),1),zeros(length(UnitMdlAccPavSucrose(:,4)),1));
+shuffd=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),zeros(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccPavAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1),ones(length(UnitMdlAccDSSucrose(:,4)),1),zeros(length(UnitMdlAccDSSucrose(:,4)),1),ones(length(UnitMdlAccPavSucrose(:,4)),1),zeros(length(UnitMdlAccPavSucrose(:,4)),1));
+rat=cat(1,RvpppasALL.Rat,RvpppasALL.Rat,RvppasALL.Rat,RvppasALL.Rat,RPPSall.Rat,RPPSall.Rat,RPSall.Rat,RPSall.Rat);
+exposure=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccDSAlcohol(:,4)),1),ones(length(UnitMdlAccPavAlcohol(:,4)),1),ones(length(UnitMdlAccPavAlcohol(:,4)),1),zeros(length(UnitMdlAccDSSucrose(:,4)),1),zeros(length(UnitMdlAccDSSucrose(:,4)),1),zeros(length(UnitMdlAccPavSucrose(:,4)),1),zeros(length(UnitMdlAccPavSucrose(:,4)),1));
+tbl=table(cat(1,UnitMdlAccDSAlcohol(:,4),UnitMdlAccDSShuffAlcohol(:,4),UnitMdlAccPavAlcohol(:,4),UnitMdlAccPavShuffAlcohol(:,4),UnitMdlAccDSSucrose(:,4),UnitMdlAccDSShuffSucrose(:,4),UnitMdlAccPavSucrose(:,4),UnitMdlAccPavShuffSucrose(:,4)),task,shuffd,rat,exposure,'variablenames',{'Accuracy','Task','Shuffled','Subject','Exposure'});
+lme=fitlme(tbl,'Accuracy~Task+Exposure+Task:Exposure+(1|Subject)');%+(1|Subject:Task:Exposure)');
+
+% [~,AccAnova,AccStats]=anovan(cat(1,UnitMdlAccDSAlcohol(:,4),UnitMdlAccDSShuffAlcohol(:,4),UnitMdlAccPavAlcohol(:,4),UnitMdlAccPavShuffAlcohol(:,4)),{task,shuffd,rat},'varnames',{'task','shuff vs true','subject'},'random',[3],'nested',[0 0 0;0 0 0;1 0 0],'model','full');
+% %no subject
+% [~,AccAnova,AccStats]=anovan(cat(1,UnitMdlAccDSAlcohol(:,4),UnitMdlAccDSShuffAlcohol(:,4),UnitMdlAccPavAlcohol(:,4),UnitMdlAccPavShuffAlcohol(:,4)),{task,shuffd},'varnames',{'task','shuff vs true',},'model','full');
+% %no shuff vs true
+% fitlme
+% task=cat(1,ones(length(UnitMdlAccDSAlcohol(:,4)),1),zeros(length(UnitMdlAccPavAlcohol(:,4)),1));
+% rat=cat(1,RvpppasALL.Rat,RvppasALL.Rat);
+% [~,AccAnova,AccStats]=anovan(cat(1,UnitMdlAccDSAlcohol(:,4),UnitMdlAccPavAlcohol(:,4)),{task,rat},'varnames',{'task','subject'},'random',[2],'nested',[0 0;1 0],'model','full');
+
+
+
+