%% % 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) & Tmpostwin(k,1) & Tmprewin(k,1))/(prewin(k,2)-prewin(k,1)); ev(k).post(m)=sum(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)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