3Tasks_paper/Scripts/A_DataProcessing_FR5.m

%% 

clear all; clc; 
global Dura Baseline Tm Tbase BSIZE Tbin DuraITI TmITI
tic
path='E:\Youna\Post-doc\MATLAB\Youna Matlab\DT Nex Files\RESULTfr5.xls'; %excel file with windows etc


%Main settings
BSIZE=0.01; %Do not change
Dura=[-25 20]; Tm=Dura(1):BSIZE:Dura(2);
DuraITI=[-50 20]; TmITI=DuraITI(1):BSIZE:DuraITI(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.01; %for comparing pre versus post windows, or event A versus event B
MinNumTrials=10;
R=[];R.Ninfo={};NN=0;
Nneurons=[];
Neur=0;
BinSize=0.05;

load('RAW_FR5.mat');
RAW=RAW;

%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:a8'); % cell that contains the event names
prewin =  [-1 0];
EventWin =  xlsread(path,'Windows','b3:c8'); %for FR5
PreEventWin = xlsread(path,'Windows','d3:e8');
PostEventWin = xlsread(path,'Windows','f3:g8');
[~,Session] = xlsread(path,'Windows','c20:c29');

%% Pre-Allocate

R.Erefnames= Erefnames;
R.Session= Session;
SessionName = {RAW.SessionName};
Ncum=0;

%Finds the total number of neurons accross all sessions
for i=1:length(Session) %loops through sessions
    Ses=strcmp(Session(i),SessionName);
    SesIndex=find(Ses==1);    
    NN=0; Neur=0;
    R.Ses(i).Ninfo=[];
    
    for y=1:length(SesIndex)
        R.Ses(i).Ninfo=cat(1,R.Ses(i).Ninfo,RAW(SesIndex(y)).Ninfo);
        Neur=Neur+size(RAW(SesIndex(y)).Nrast,1);
        Nneurons(i)=Neur;
    end
    
        load('Coord_FR5.mat');
        R.Ses(i).Coord=part(i).Coord;
        R.Ses(i).Structure=part(i).Coord(:,4);
        R.Ses(i).Ninfo=cat(2, R.Ses(i).Ninfo, mat2cell([1:Nneurons(i)]',ones(1,Nneurons(i)),1));
        
    for k=1:length(Erefnames)
        R.Ses(i).Ev(k).PSTHraw(1:Nneurons(i),1:length(Tm))=NaN(Nneurons(i),length(Tm));
        R.Ses(i).Ev(k).PSTHz(1:Nneurons(i),1:length(Tm))=NaN(Nneurons(i),length(Tm));
        R.Ses(i).Ev(k).PSTHrawBL(1:Nneurons(i),1:length(Tbase))=NaN(Nneurons(i),length(Tbase));
        R.Ses(i).Ev(k).Meanraw(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).rawMeanz(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).rawMeanzPre(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).Meanz(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).ttestPreEvent(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).ttestPostEvent(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).RespDirPre(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).RespDirPost(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).NumberTrials(1:Nneurons(i),1)=NaN;
    end   
end


% preallocating
R.Param.Tm=Tm;
R.Param.Tbin=Tbin;
R.Param.Dura=Dura;
R.Param.Baseline=Baseline;
R.Param.PStat=PStat;
R.Param.MinNumTrials=MinNumTrials;
R.Param.path=path;
R.Param.prewin=prewin;
R.Param.postwinPreEvent=PreEventWin;
R.Param.postwinPostEvent=PostEventWin;
R.Param.SmoothTYPE=SmoothTYPE;
R.Param.SmoothSPAN=SmoothSPAN;

%% runs the main routine

NS=0;
TNN=0;              
for i=1:length(Session) %loops through sessions
    Ses=strcmp(Session(i),SessionName);
    SesIndex=find(Ses==1);    
    NN=0;
    NS=NS+1;
    
    for y=1:length(SesIndex)
       for j= 1:size(RAW(SesIndex(y)).Nrast,1) %Number of neurons per session            
            NN=NN+1; %neuron counter
            TNN=TNN+1; %neuron counter across all session
 
            if R.Ses(i).Structure(NN)~=0 %to avoid analyzing excluded neurons
                
                ev2=length(RAW(SesIndex(y)).Nrast{j})/RAW(SesIndex(y)).Nrast{j}(end);% FR computed across whole session.
                R.WF(TNN,:)=RAW(SesIndex(y)).waveforms(j,:);
                R.WFanalysis(TNN,1)=ev2; 
                proportionISI=0;
                for k=3:1:size(RAW(SesIndex(y)).Nrast{j,1},1) % loop thru timestamps waveform
                    ISIn=RAW(SesIndex(y)).Nrast{j,1}(k)-RAW(SesIndex(y)).Nrast{j,1}(k-1); 
                    ISIo=RAW(SesIndex(y)).Nrast{j,1}(k-1)-RAW(SesIndex(y)).Nrast{j,1}(k-2);
                    calcul(k-2)=2*abs(ISIn-ISIo)/(ISIn+ISIo);
                    if ISIo>0.5
                        proportionISI=proportionISI+1;
                    end
                end
                R.CV2(TNN,1)=mean(calcul);%CV2
                R.CV2(TNN,2)=100*proportionISI/(size(RAW(SesIndex(y)).Nrast{j,1},1)-1);% of ISI>0.5s
                R.SESSION(TNN,1)=NS;
                
                LP1Ind=strcmp('First_LP',RAW(SesIndex(y)).Einfo(:,2)); %find the event id number from RAW
                for k=1:length(Erefnames) %loops thorough the events
                    EvInd=strcmp(Erefnames(k),RAW(SesIndex(y)).Einfo(:,2)); %find the event id number from RAW
                    if sum(EvInd)==0
                        fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
                    end
                    if k==1
                        FirstPressInd=EvInd;                     
                    end
                    
                    
                    R.Ses(i).Ev(k).NumberTrials(NN,1)=length(RAW(SesIndex(y)).Erast{EvInd});
                    Tbase=prewin(1,1):BSIZE:prewin(1,2); % for FR5
                    [PSR0,N0]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{LP1Ind},prewin(1,:),{1});
                    [PSR1,N1]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)                    
                    [PSR2,N2]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
                    BLcell=MakePSR04(RAW(SesIndex(y)).Erast(FirstPressInd),RAW(SesIndex(y)).Erast{EvInd},[-1800 1],{2,'Last'});
                    
                    if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
                        %Fixed bin size
                    [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(SesIndex(y)).Nrast{j},1),1),{2,0,BinSize});%-----DP used here
                    [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(SesIndex(y)).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
                        % calculate MeanZ from PSTH before smoothing
                        
                        if sum(PTH0,2)~=0
                            unsmoothPSTHz(1,1:length(Tm))=normalize(PTH1,PTH0,0);
                            R.Ses(i).Ev(k).rawMeanz(NN,1)=nanmean(unsmoothPSTHz(1,Tm>EventWin(k,1) & Tm<EventWin(k,2)),2);
                            R.Ses(i).Ev(k).rawMeanzPre(NN,1)=nanmean(unsmoothPSTHz(1,Tm>PreEventWin(k,1) & Tm<PreEventWin(k,2)),2);
                        else
                            R.Ses(i).Ev(k).rawMeanz(NN,1)=NaN;
                            R.Ses(i).Ev(k).rawMeanzPre(NN,1)=NaN;
                        end
                            
                        PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
                        PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';
                        
                        %------------- Fills the R.Ses(i).Ev(k) fields --------------
                        R.Ses(i).Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
                        R.Ses(i).Ev(k).PSTHrawBL(NN,1:length(Tbase))=PTH0;
                        R.Ses(i).Ev(k).Meanraw(NN,1)=nanmean(R.Ses(i).Ev(k).PSTHraw(NN,Tm>PostEventWin(k,1) & Tm<PostEventWin(k,2)),2);
                        if sum(PTH0,2)~=0
                            R.Ses(i).Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
                            R.Ses(i).Ev(k).Meanz(NN,1)=nanmean(R.Ses(i).Ev(k).PSTHz(NN,Tm>EventWin(k,1) & Tm<EventWin(k,2)),2);
                        else
                            R.Ses(i).Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
                            R.Ses(i).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(SesIndex(y)).Erast{EvInd},1),1);
                        ev(k).PreEvent=NaN(size(RAW(SesIndex(y)).Erast{EvInd},1),1);
                        ev(k).PostEvent=NaN(size(RAW(SesIndex(y)).Erast{EvInd},1),1);
                        for m=1:size(RAW(SesIndex(y)).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).pre(m)=sum(PSR2{m}<(BLcell{1,m}+prewin(1,2)) & PSR2{m}>(BLcell{1,m}+prewin(1,1)))/(prewin(1,2)-prewin(1,1));
                           ev(k).PreEvent(m)=sum(PSR2{m}<PreEventWin(k,2) & PSR2{m}>PreEventWin(k,1))/(PreEventWin(k,2)-PreEventWin(k,1));
                           ev(k).PostEvent(m)=sum(PSR2{m}<PostEventWin(k,2) & PSR2{m}>PostEventWin(k,1))/(PostEventWin(k,2)-PostEventWin(k,1));
                        end

                        %---------------------------Response detection w/ SignRank on pre/post windows---------------------------
                        
                        if  ~isempty(EvInd) && R.Ses(i).Ev(k).NumberTrials(NN,1)>= MinNumTrials  && nanmean(ev(k).pre(:,1),1)>0.5 %avoid analyzing sessions where that do not have enough trials%avoid analyzing sessions where that do not have enough trials
                            
                            [~,R.Ses(i).Ev(k).ttestPreEvent(NN,1)]=signrank(ev(k).pre, ev(k).PreEvent); %Signrank used here because it is a dependant sample test
                            if R.Ses(i).Ev(k).ttestPreEvent(NN,1)<PStat
                                R.Ses(i).Ev(k).RespDirPre(NN,1)=sign(mean(ev(k).PreEvent)-mean(ev(k).pre));                               
                            else R.Ses(i).Ev(k).RespDirPre(NN,1)=0;
                            end
                            
                            [~,R.Ses(i).Ev(k).ttestPostEvent(NN,1)]=signrank(ev(k).pre, ev(k).PostEvent); %Signrank used here because it is a dependant sample test
                            if R.Ses(i).Ev(k).ttestPostEvent(NN,1)<PStat
                                R.Ses(i).Ev(k).RespDirPost(NN,1)=sign(mean(ev(k).PostEvent)-mean(ev(k).pre));                          
                            else R.Ses(i).Ev(k).RespDirPost(NN,1)=0;
                            end
                            
                            if R.Ses(i).Ev(k).RespDirPre(NN,1)>0 || R.Ses(i).Ev(k).RespDirPost(NN,1) >0
                                Search=find(Tm>=PreEventWin(k,1) & Tm<=PostEventWin(k,2));
                                [R.Ses(i).Ev(k).MaxVal(NN,1),MaxInd]=max(R.Ses(i).Ev(k).PSTHraw(NN,Search));
                                R.Ses(i).Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                            else
                                Search=find(Tm>=PreEventWin(k,1) & Tm<PostEventWin(k,2));
                                [R.Ses(i).Ev(k).MaxVal(NN,1),MaxInd]=min(R.Ses(i).Ev(k).PSTHraw(NN,Search));
                                R.Ses(i).Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                            end
                        else R.Ses(i).Ev(k).RespDirEvent(NN,1)=0;
                        end
                            
                            
                        end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
                    end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
                %Events
                
                
                fprintf('Neuron ID # %d\n',NN);
               elseif R.Ses(i).Structure(NN)~=0
            end %exclusion: IF R.Structure(NN)~=0 to avoid analyzing excluded neurons
        end %neurons: FOR j= 1:size(RAW(SesIndex(y)).Nrast,1)
    end %SesIndex
end %sessions: FOR i=1:length(RAW)

for k=1:size(R.Ses,2)
    tableTRN=[];
    for i=1:length(Erefnames)
        tableTRN(:,i)=R.Ses(k).Ev(i).RespDirPre(:,1);
        tableTRN(:,i+length(Erefnames))=R.Ses(k).Ev(i).RespDirPost(:,1);
    end
    R.Ses(k).TRN(:,1)=nansum(abs(tableTRN),2);
end

save('R_FR5.mat','R')

toc