123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213 |
- function [DATout]=ResDetect08(Rin,Rbasein,INTin,NORMTYPin,WINin,CI)
- % This function detects responses in peri-event histograms. Confidence intervals are extracted from the
- % distribution of the baseline period.
- % Arpil 2012, Fred: This version was adapted from RespDetect07 (now handles single-neuron PTH)
- % This version contains the helper functions "sigtestin", "normalizein", "findonset03in"
- % Inputs:
- % Rin is a matrix e.g. PTH1(neurons, bins)
- % Rbasein is the baseline version of Rin (used to normalize)
- % INTin is the window in which responses are searched
- % NORMTYPin for the function normalizein - should be set to zero
- % WINin are onset and offset requirements; by default it is [0.1 0.3]
- % CI is the confidence interval (inhibitions<lower bound, excitation>Upper bound)
- % Outputs:
- % DATout.R= Rin
- % .LUb(neurons, 1:2): Lower and Upper Bounds = Threshold ofr excitation and inhibitions
- % .Rnorm(neurons, bins)= z-score transformed Rin
- % .Rnormmask(neurons, bins) = z-score with all bins that are within the confidence intervals are zeroed
- % .RFLAG(neurons, 1): +1 for excitation, -1 for inhibition, 0 for no resp.
- % .CFLAG(neurons, 1)= 0 or 1. 1 means that an inhibition follows an excation or vice-versa
- % .Onsets(neurons,1:2): First and Second columns stores the onset of excitations and inhibitions respectively
- % .Offsets(neurons,1:2): same as onset
- % .Paramnames={'Dura','Baseline','Bsize','INT','WINs'};
- % .Params={Dura,Baseline,BSIZE,INTin,WINin};
- global Dura Tm Tbase BSIZE Erefnames
- %Allows to handle single neurones cases where R(bins, 1) instead of R(1,neurons)
- if size(Rin,2)==1
- R=Rin'; Rbase=Rbasein';
- else
- R=Rin;Rbase=Rbasein;
- end
- if isempty(WINin)
- WINON=ceil(0.1/BSIZE);WINOFF=ceil(0.3/BSIZE);
- WINFLAG=0;
- elseif size(WINin,1)==1
- WINON=ceil(WINin(1)/BSIZE);WINOFF=ceil(WINin(2)/BSIZE); %same detection for all neurons
- WINFLAG=0;
- else
- WINONmat=ceil(WINin(:,1)/BSIZE);WINOFFmat=ceil(WINin(:,2)/BSIZE); % different for each neuron
- WINFLAG=1;
- end
- if size(INTin,1)==1
- Iref=find(Tm>=INTin(1) & Tm<=INTin(2));
- INTFLAG=0;
- else
- INTmat=INTin;
- INTFLAG=1;
- end
- % ------------ Determines the thresholds and makes a mask ------------
- %ALPHA=0.05/(size(R,2)-length(Tmind)); %Bonferroni correction
- if length(CI)==1, CI0(1)=CI;CI0(2)=CI; else CI0=CI; end %size(CI0)=(1x2)
- if CI0(1)>0
- %LUb= Lower and Upper Bound = threshold determined by the confidence intervals
- LUb=(prctile(Rbase',[CI0(1)/2 100-CI0(2)/2]))'; % Outside LUb range consider significant
- if size(LUb,2)==1, LUb=LUb'; end %To handle single neuron cases
- else
- %Ali's stuff - I think it uses it to go beyond 0 & 100% confidence intervals
- %it can be reached by using negative CI values
- LUb(:,1)=min(Rbase,[],2)-(max(Rbase,[],2)-min(Rbase,[],2))*abs(CI0(1))/200;
- LUb(:,2)=max(Rbase,[],2)+(max(Rbase,[],2)-min(Rbase,[],2))*abs(CI0(2))/200;
- end
- [lmask,umask]=sigtestin(R,LUb);
- Rnorm=normalizein(R,Rbase,NORMTYPin);%used to put (-) to bins below lower bound
- Rnormmask=Rnorm.*(umask+lmask);
- %Outupts R, LUb, Normmask and Rnorm
- DATout.R=R;DATout.LUb=LUb;DATout.Rnormmask=Rnormmask;DATout.Rnorm=Rnorm;
- for j=1:size(R,1) %loops through neurons
- if mod(j,11)==0
- fprintf('.');
- elseif mod(j,100)==0
- fprintf('%d',j)
- end
- if WINFLAG
- WINON=WINONmat(j);WINOFF=WINOFFmat(j);
- end
- if INTFLAG
- Iref=find(Tm>INTmat(j,1) & Tm<INTmat(j,2));
- end
- %PInd1 and NInd1 are the Tm indexes of the first excitation and inhibition respectively
- PInd1=findonset03in(double(Rnormmask(j,Iref)>0),1,WINON);
- NInd1=findonset03in(double(Rnormmask(j,Iref)<0),1,WINON);
-
-
- %% Flag settings
- %------------ Finds the onsets/offset of the FIRST response detected ------------
- Ponset=NaN;Nonset=NaN;Poffset=NaN;Noffset=NaN;TFLAG=[0 0];%termination flag
- if PInd1<NInd1
- % RFLAG= Response Flag: +1:excitation, -1:inhibition, 0:No response
- RFLAG=1;
- % Ponset is the excitation onset
- Ponset=Tm(Iref(PInd1));
- elseif PInd1>NInd1
- RFLAG=-1;
- % Nonset is the inhibition onset
- Nonset=Tm(Iref(NInd1));
- elseif isinf(PInd1) && isinf(NInd1)
- RFLAG=0;
- end
-
- % CFLAG= conflict flag=1 if both excitation and inhibition detected
- if ~isinf(PInd1) && ~isinf(NInd1)
- CFLAG=1;
- Ponset=Tm(Iref(PInd1));Nonset=Tm(Iref(NInd1));
- else
- CFLAG=0;
- end
-
- %------------ Finds the onsets/offset of the SECOND response detected ------------
- if ~isinf(PInd1)
- PInd2=findonset03in(double(Rnormmask(j,:)<=0),Iref(PInd1)+1,WINOFF);
- if ~isinf(PInd2)
- Poffset=Tm(PInd2);
- if Poffset<-0.03 % 0.025 is the interp1 offset with 0.05 bin
- TFLAG(1)=1;
- end
- else
- Poffset=Tm(end)+sqrt(-1); % if offset not detected it is beyond Tm(end)
- PInd2=length(Tm);
- end
- end
- if ~isinf(NInd1)
- NInd2=findonset03in(double(Rnormmask(j,:)>=0),Iref(NInd1)+1,WINOFF);
- if ~isinf(NInd2)
- Noffset=Tm(NInd2);
- if Noffset<-0.03 % 0.025 is the interp1 offset with 0.05 bin
- TFLAG(2)=1;
- end
- else
- Noffset=Tm(end)+sqrt(-1); % if offset not detected it is beyond Tm(end)
- NInd2=length(Tm);
- end
- end
-
- % setting outputs (new change to work for LP,Ent and Ext while keeping DS,NS results intact
- if (CFLAG==0 && sum(TFLAG)==1) || (CFLAG==1 && sum(TFLAG)==2) % early termination response not valid (set as imaginary)
- DATout.RFLAG(j,1)=RFLAG*sqrt(-1);
- else
- DATout.RFLAG(j,1)=RFLAG;
- end
- DATout.RFLAG(j,1)=RFLAG;
- DATout.CFLAG(j,1)=CFLAG;
- % offset and onsets should be swaped and negated
- DATout.Onsets(j,:)=[Ponset,Nonset];DATout.Offsets(j,:)=[Poffset,Noffset];
- end
- % setting outputs
- DATout.Paramnames={'Dura','Bsize','INT','WINs'};
- DATout.Params={Dura,BSIZE,INTin,WINin};
- %fprintf('\n response detection completed \n')
- %----- HELPER FUNCTIONS ----
- function [LMask,UMask]=sigtestin(Matin,Bin)
- UMask=zeros(size(Matin));LMask=zeros(size(Matin));
- for k=1:size(Matin,1)
- TMP=Matin(k,:);B=Bin(k,:);
- LMask(k,find(TMP<(B(1)-eps)))=1;
- UMask(k,find(TMP>(B(2)+eps)))=1;
- end
- function Mnormal=normalizein(M,Mbase,MD)
- %MD is the type of normalization.
- % 0 = Z-score >>(X-mean)/SD
- % 1 = normalized to response max
- % 2 = Baseline substracted
- % 3 = normalized by SD only >> X/SD
- %Added April 2012 to handle vectors
- if size(M,2)==1, M=M'; Mbase=Mbase'; end
- Mnormal=zeros(size(M));
- for k=1:size(M,1)
- if MD==0
- Mnormal(k,:)=(M(k,:)-nanmean(Mbase(k,:)))/nanstd(Mbase(k,:)); % Zscore
- elseif MD==1
- Tmp=M(k,:)-nanmean(Mbase(k,:));
- Mnormal(k,:)=Tmp/max(abs(Tmp));% mean corrected normalized to RESPONSE MAX
- elseif MD==2
- Mnormal(k,:)=(M(k,:)-nanmean(Mbase(k,:)));
- elseif MD==3
- Mnormal(k,:)=M(k,:)/nanstd(Mbase(k,:)); % Zscore no demean
- end
- end
- function Indout=findonset03in(Rin,Indin,Dur)
- % Finds the index of the first element when 'Dur' consecutive elements are equal to 1
- % Rin must contain double and not logical values. Therefore, the result of a find function
- % must be converted as follow: double(find(MAT>X))
- % Indin: looks from index Indin onward
- % Dur: Indetfies onset only if response is Dur long
- Indout=Inf;
- if ~isinf(Indin)
- R=Rin(Indin:end);
- TMP=find(conv(ones(1,Dur),R)==Dur)-(Dur-1);
- if ~isempty(TMP)
- Indout=TMP(1)+Indin-1;
- end
- end
|