123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333 |
- %%
- % 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
|