VentralPallidumAlcohol2019/AnalysisAndFigures/e_SingleUnitPavlovianAlcoholAnalysis.m

%%
% 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='RESULTpa.xls';
load('RAWvppaALL.mat')
RAW=RAWvppaALL;

%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 - this code uses a fixed bin width
PStat=0.05; %for comparing pre versus post windows, or event A versus event B
MinNumTrials=3;
R=[];RvppaALL.Ninfo={};NN=0;Nneurons=0;
BinSize=0.02;

%Smoothing
Smoothing=0; %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)
    RvppaALL.Ninfo=cat(1,RvppaALL.Ninfo,RAW(i).Ninfo);
    Nneurons=Nneurons+size(RAW(i).Nrast,1);
end
load('CoordPAall.mat');
RvppaALL.Coord=CoordPAall;
%RvppaALL.Coord=ones(Nneurons,4); %comment this line and uncomment the 2 lines above when you have real coordinates.
RvppaALL.Structure=RvppaALL.Coord(:,4);
RvppaALL.Ninfo=cat(2, RvppaALL.Ninfo, mat2cell([1:Nneurons]',ones(1,Nneurons),1));
RvppaALL.Erefnames= Erefnames;

RvppaALL.Rat(1:Nneurons,1)=NaN;
RvppaALL.Session(1:Nneurons,1)=NaN;

for i=1:length(RvppaALL.Ninfo)
    RvppaALL.Rat(i,1)=str2num(RvppaALL.Ninfo{i,1}(3:5));
    RvppaALL.Session(i,1)=str2num(RvppaALL.Ninfo{i,1}(10:11));
end

% preallocating
RvppaALL.Param.Tm=Tm;
RvppaALL.Param.Tbin=Tbin;
RvppaALL.Param.Dura=Dura;
RvppaALL.Param.Baseline=Baseline;
RvppaALL.Param.PStat=PStat;
RvppaALL.Param.MinNumTrials=MinNumTrials;
RvppaALL.Param.path=path;
RvppaALL.Param.prewin=prewin;
RvppaALL.Param.postwin=postwin;
RvppaALL.Param.BLwin=BLWin;
RvppaALL.Param.RespWin=RespWin;
RvppaALL.Param.ResponseReq=WINin;
RvppaALL.Param.SmoothTYPE=SmoothTYPE;
RvppaALL.Param.SmoothSPAN=SmoothSPAN;


for k=1:length(Erefnames)
    RvppaALL.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    RvppaALL.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    RvppaALL.Ev(k).Meanraw(1:Nneurons,1)=NaN;
    RvppaALL.Ev(k).Meanz(1:Nneurons,1)=NaN;
    RvppaALL.Ev(k).BW(1:Nneurons,1)=NaN;
    RvppaALL.Ev(k).ttest(1:Nneurons,1)=NaN;
    RvppaALL.Ev(k).RespDir(1:Nneurons,1)=NaN;
    RvppaALL.Ev(k).RFLAG(1:Nneurons,1)=NaN;
    RvppaALL.Ev(k).CFLAG(1:Nneurons,1)=NaN;
    RvppaALL.Ev(k).Onsets(1:Nneurons,:)=NaN(Nneurons,2);
    RvppaALL.Ev(k).Offsets(1:Nneurons,:)=NaN(Nneurons,2);
    RvppaALL.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
    RvppaALL.Ev(k).MaxVal(1:Nneurons,1)=NaN;
    RvppaALL.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 RvppaALL.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
                
                RvppaALL.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
                Tbase=prewin(k,1):BSIZE:prewin(k,2);
                if  ~isempty(EvInd) && RvppaALL.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 RvppaALL.Ev(k) fields --------------
                        RvppaALL.Ev(k).BW(NN,1)=BW1;
                        RvppaALL.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
                        RvppaALL.Ev(k).Meanraw(NN,1)=nanmean(RvppaALL.Ev(k).PSTHraw(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
                        if sum(PTH0,2)~=0
                            RvppaALL.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
                            RvppaALL.Ev(k).Meanz(NN,1)=nanmean(RvppaALL.Ev(k).PSTHz(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
                        else
                            RvppaALL.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
                            RvppaALL.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---------------------------
                        [~,RvppaALL.Ev(k).ttest(NN,1)]=ttest(ev(k).pre, ev(k).post); %Signrank used here because it is a dependant sample test
                        if RvppaALL.Ev(k).ttest(NN,1)<PStat
                            RvppaALL.Ev(k).RespDir(NN,1)=sign(mean(ev(k).post)-mean(ev(k).pre));
                            if RvppaALL.Ev(k).RespDir(NN,1)>0
                                Search=find(Tm>=postwin(k,1) & Tm<=postwin(k,2));
                                [RvppaALL.Ev(k).MaxVal(NN,1),MaxInd]=max(RvppaALL.Ev(k).PSTHraw(NN,Search));
                                RvppaALL.Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                            else
                                Search=find(Tm>=postwin(k,1) & Tm<=postwin(k,2));
                                [RvppaALL.Ev(k).MaxVal(NN,1),MaxInd]=min(RvppaALL.Ev(k).PSTHraw(NN,Search));
                                RvppaALL.Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                            end
                        else RvppaALL.Ev(k).RespDir(NN,1)=0;
                        end
                        
                        %---------------------------Response detection w/ RESDETECT08 on running windows---------------------------
                        RD=ResDetect08(PTH1,PTH0,RespWin(k,:),0,WINin, CI);
                        
                        RvppaALL.Ev(k).RFLAG(NN,1)=RD.RFLAG;
                        RvppaALL.Ev(k).CFLAG(NN,1)=RD.CFLAG;
                        RvppaALL.Ev(k).Onsets(NN,:)=RD.Onsets';
                        RvppaALL.Ev(k).Offsets(NN,:)=RD.Offsets';
                        RvppaALL.Param.Paramnames=RD.Paramnames;
                        RvppaALL.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(CSPlus)~=0 & sum(CSMinus)~=0
                    RvppaALL.CScriteria=linspace(0, max([ev(CSPlus).post(:);ev(CSMinus).post(:)]),200);%Determine the criteria range
                    RvppaALL.CSprecriteria=linspace(0, max([ev(CSPlus).pre(:);ev(CSMinus).pre(:)]),200);%Determine the criteria range
                    for s=1:length(RvppaALL.CScriteria);
                        RvppaALL.CSfalsepos(NN,s)=sum(ev(CSMinus).post>=RvppaALL.CScriteria(s))/length(ev(CSMinus).post); %Determine the probability that the baseline is greater than each criteria (x)
                        RvppaALL.CStruepos(NN,s)=sum(ev(CSPlus).post>=RvppaALL.CScriteria(s))/length(ev(CSPlus).post);%Determine the probability that the post-window firing firing is greater than critera(y)
                        RvppaALL.CSprefalsepos(NN,s)=sum(ev(CSMinus).pre>=RvppaALL.CScriteria(s))/length(ev(CSMinus).pre); %Determine the probability that the baseline is greater than each criteria (x)
                        RvppaALL.CSpretruepos(NN,s)=sum(ev(CSPlus).pre>=RvppaALL.CScriteria(s))/length(ev(CSPlus).pre);%Determine the probability that the pre-window firing firing is greater than critera(y)
                    end
                    RvppaALL.CSauROC(NN,1)=trapz(RvppaALL.CSfalsepos(NN,:),RvppaALL.CStruepos(NN,:));%Calculate area under the curve for this neuron for this event
                    RvppaALL.CSpreauROC(NN,1)=trapz(RvppaALL.CSprefalsepos(NN,:),RvppaALL.CSpretruepos(NN,:));%Calculate area under the curve for this neuron for this event
                    [RvppaALL.CSStat(NN,1),~]=ranksum(ev(CSPlus).post,ev(CSMinus).post); %Ranksum test used becasue it is an independant sample test
                else RvppaALL.CSStat(NN,1)=NaN;
                    RvppaALL.CSauROC(NN,1)=NaN;
                    RvppaALL.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 RvppaALL.Ev(PECSPlus).RespDir(NN,1)~=0 || RvppaALL.Ev(PECSMinus).RespDir(NN,1)~=0
            %                     [RvppaALL.PECueStat(NN,1),~]=ranksum(ev(PECSPlus).post,ev(PECSMinus).post); %Ranksum test used becasue it is an independant sample test
            %                 else RvppaALL.PECueStat(NN,1)=NaN;
            %                 end
            %             end
            
            FirstPE=strcmp('FirstPE', Erefnames);
            PortEntry=strcmp('PortEntry', Erefnames);
            if sum(FirstPE)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
                if RvppaALL.Ev(FirstPE).RespDir(NN,1)~=0 || RvppaALL.Ev(PortEntry).RespDir(NN,1)~=0
                    [RvppaALL.PERewStat(NN,1),~]=ranksum(ev(FirstPE).post,ev(PortEntry).post); %Ranksum test used becasue it is an independant sample test
                else RvppaALL.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 RvppaALL.Ev(PECSMinus).RespDir(NN,1)~=0 || RvppaALL.Ev(ITIPE).RespDir(NN,1)~=0
            %                     [RvppaALL.PEUnrewCueStat(NN,1),~]=ranksum(ev(PECSMinus).post,ev(ITIPE).post); %Ranksum test used becasue it is an independant sample test
            %                 else RvppaALL.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(CSPlusnoPE)~=0 & sum(CSPluswPE)~=0
                    RvppaALL.CSrespcriteria=linspace(0, max([ev(CSPlusnoPE).post(:);ev(CSPluswPE).post(:)]),200);%Determine the criteria range
                    RvppaALL.CSprerespcriteria=linspace(0, max([ev(CSPlusnoPE).pre(:);ev(CSPluswPE).pre(:)]),200);%Determine the criteria range
                    for s=1:length(RvppaALL.CSrespcriteria);
                        RvppaALL.CSrespfalsepos(NN,s)=sum(ev(CSPlusnoPE).post>=RvppaALL.CSrespcriteria(s))/length(ev(CSPlusnoPE).post); %Determine the probability that the baseline is greater than each criteria (x)
                        RvppaALL.CSresptruepos(NN,s)=sum(ev(CSPluswPE).post>=RvppaALL.CSrespcriteria(s))/length(ev(CSPluswPE).post);%Determine the probability that the post-window firing firing is greater than critera(y)
                        RvppaALL.CSprerespfalsepos(NN,s)=sum(ev(CSPlusnoPE).pre>=RvppaALL.CSrespcriteria(s))/length(ev(CSPlusnoPE).pre); %Determine the probability that the baseline is greater than each criteria (x)
                        RvppaALL.CSpreresptruepos(NN,s)=sum(ev(CSPluswPE).pre>=RvppaALL.CSrespcriteria(s))/length(ev(CSPluswPE).pre);%Determine the probability that the post-window firing firing is greater than critera(y)
                    end
                    RvppaALL.CSrespauROC(NN,1)=trapz(RvppaALL.CSrespfalsepos(NN,:),RvppaALL.CSresptruepos(NN,:));%Calculate area under the curve for this neuron for this event
                    RvppaALL.CSprerespauROC(NN,1)=trapz(RvppaALL.CSprerespfalsepos(NN,:),RvppaALL.CSpreresptruepos(NN,:));
                    [RvppaALL.CSPlusResponseStat(NN,1),~]=ranksum(ev(CSPluswPE).post,ev(CSPlusnoPE).post); %Ranksum test used becasue it is an independant sample test
                else
                    RvppaALL.CSrespcriteria=linspace(NaN,NaN,200);
                    RvppaALL.CSrespauROC(NN,1)=NaN;
                    RvppaALL.CSPlusResponseStat(NN,1)=NaN;
                    RvppaALL.CSprerespauROC(NN,1)=NaN;
                    RvppaALL.CSrespfalsepos(NN,1:length(RvppaALL.CSrespcriteria))=NaN;
                    RvppaALL.CSrespfalsepos(NN,1:length(RvppaALL.CSrespcriteria))=NaN;
                    fprintf('NAN! Neuron ID # %d\n',NN);
                end
            else
                RvppaALL.CSrespcriteria=linspace(NaN,NaN,200);
                RvppaALL.CSrespauROC(NN,1)=NaN;
                RvppaALL.CSPlusResponseStat(NN,1)=NaN;
                RvppaALL.CSprerespauROC(NN,1)=NaN;
                RvppaALL.CSrespfalsepos(NN,1:length(RvppaALL.CSrespcriteria))=NaN
                RvppaALL.CSrespfalsepos(NN,1:length(RvppaALL.CSrespcriteria))=NaN;
                fprintf('NAN! Neuron ID # %d\n',NN);
            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 RvppaALL.Ev(CSMinuswPE).RespDir(NN,1)~=0 || RvppaALL.Ev(CSMinusnoPE).RespDir(NN,1)~=0
            %                     [RvppaALL.CSMinusResponseStat(NN,1),~]=ranksum(ev(CSMinuswPE).post,ev(CSMinusnoPE).post); %Ranksum test used becasue it is an independant sample test
            %                 else RvppaALL.CSMinusResponseStat(NN,1)=NaN;
            %                 end
            %             end
            %
          
    RvppaALL.CSPlusRatio(NN,1)=RvppaALL.Ev(3).NumberTrials(NN)/length(ev(CSPlus).post);
    RvppaALL.CSMinusRatio(NN,1)=RvppaALL.Ev(5).NumberTrials(NN)/length(ev(CSMinus).post);
            
            fprintf('Neuron ID # %d\n',NN);
        %elseif RvppaALL.Structure(NN)~=0
        %end %exclusion: IF RvppaALL.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)

RvppaALL.CSfalseposMEAN(1,:)=nanmean(RvppaALL.CSfalsepos(RvppaALL.Structure==10,:));
RvppaALL.CSfalseposSEM(1,:)=nanste(RvppaALL.CSfalsepos(RvppaALL.Structure==10,:),1);

RvppaALL.CStrueposMEAN(1,:)=nanmean(RvppaALL.CStruepos(RvppaALL.Structure==10,:));
RvppaALL.CStrueposSEM(1,:)=nanste(RvppaALL.CStruepos(RvppaALL.Structure==10,:),1);

RvppaALL.CSrespfalseposMEAN(1,:)=nanmean(RvppaALL.CSrespfalsepos(RvppaALL.Structure==10,:));
RvppaALL.CSrespfalseposSEM(1,:)=nanste(RvppaALL.CSrespfalsepos(RvppaALL.Structure==10,:),1);

RvppaALL.CSresptrueposMEAN(1,:)=nanmean(RvppaALL.CSresptruepos(RvppaALL.Structure==10,:));
RvppaALL.CSresptrueposSEM(1,:)=nanste(RvppaALL.CSresptruepos(RvppaALL.Structure==10,:),1);

if SAVE_FLAG
    save('RvppaALL.mat','RvppaALL')
end

toc