VentralPallidumAlcohol2019/AnalysisAndFigures/d_SingleUnitInstrumentalAlcoholAnalysis.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='RESULTppa.xls';
load('RAWvpppaALL.mat')
RAW=RAWvpppaALL;

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

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

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

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


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


    RvpppaALL.CSPlusRatio(NN,1)=RvpppaALL.Ev(3).NumberTrials(NN)/length(ev(CSPlus).post);
    RvpppaALL.CSMinusRatio(NN,1)=RvpppaALL.Ev(5).NumberTrials(NN)/length(ev(CSMinus).post);
    
%   
% if NN==673
%     fprintf('Neuron ID # %d\n',NN);
%     fprintf('CS Minus Response Number ')% length(ev(CSMinuswPE).post))
%     fprintf('CS Minus Number ')% length(ev(CSMinus.post)))
%     return
% else
    fprintf('Neuron ID # %d\n',NN);
% end

        %elseif RvpppaALL.Structure(NN)~=0
        %end %exclusion: IF RvpppaALL.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)

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

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

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

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

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

toc