DataAvailability_eLIFE2019/scripts/A_eventAnalysis_ExtTrain.m

%% initial analysis of event-evoked activity - Extended training dataset - neurons pooled together

% need RAW datafiles and Coord datafile with electrode placement
% need to save excel sheets in the specified folder
% the excel sheets say which events to make PSTHs for and what windows to
% use for statistical testing


%clear all; clc;

clear all; clc; 
global Dura Baseline Tm Tbase BSIZE Tbin DuraITI TmITI
tic
path='C:\Users\yvandae1\Documents\MATLAB\DT Nex Files\RESULTdt.xls'; %excel file with windows etc
load('RAWextendedtraining.mat')
RAW=RAW;

%Main settings
SAVE_FLAG=1;
FLAG_BL=0; %1 to use ITI for BL, 0 to use the pre-event window
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.05; %for comparing pre versus post windows, or event A versus event B
MinNumTrials=10;
R=[];R.Ninfo={};NN=0;Nneurons=0;
BinSize=0.050;

%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:a15'); % cell that contains the event names
prewin =  xlsread(path,'Windows','b3:c15');
PreEventWin = xlsread(path,'Windows','d3:e15');
PostEventWin = xlsread(path,'Windows','f3:g15');
EventWin =  xlsread(path,'Windows','l3:m15');
BLWin=    xlsread(path,'Windows','h3:i15');
RespWin=  xlsread(path,'Windows','j3:k15');


%Settings for Response detection w/ signrank
WINin=[0.03,0.1]; %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=[.95,.95];

%% 
%Finds the total number of neurons accross all sessions
for i=1:length(RAW)
    R.Ninfo=cat(1,R.Ninfo,RAW(i).Ninfo);
    Nneurons=Nneurons+size(RAW(i).Nrast,1);
end
load('Coord_extendedtraining.mat');
R.Coord=Coord;
%R.Coord=ones(Nneurons,4); %comment this line and uncomment the 2 lines above when you have real coordinates.
R.Structure=R.Coord(:,4); 
R.Ninfo=cat(2, R.Ninfo, mat2cell([1:Nneurons]',ones(1,Nneurons),1));
R.Sinfo=[];

for i=1:Nneurons
    fname=char(R.Ninfo(i,1));
    R.Sinfo=cat(1,R.Sinfo,cellstr(fname(11:end-4)));
end


R.Erefnames= Erefnames;



% 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.BLwin=BLWin;
R.Param.RespWin=RespWin;
R.Param.ResponseReq=WINin;
R.Param.SmoothTYPE=SmoothTYPE;
R.Param.SmoothSPAN=SmoothSPAN;
R.TRNevent(1:Nneurons,1:length(Erefnames))=NaN;


for k=1:length(Erefnames)
    R.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    R.Ev(k).PSTHrawBL(1:Nneurons,1:length(Tbase))=NaN(Nneurons,length(Tbase));
    R.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    R.Ev(k).Meanraw(1:Nneurons,1)=NaN;
    R.Ev(k).rawMeanz(1:Nneurons,1)=NaN;
    R.Ev(k).rawMeanzPre(1:Nneurons,1)=NaN;
    R.Ev(k).Meanz(1:Nneurons,1)=NaN;
    R.Ev(k).MeanzPRE(1:Nneurons,1)=NaN;
    R.Ev(k).ttestPreEvent(1:Nneurons,1)=NaN;
    R.Ev(k).ttestPostEvent(1:Nneurons,1)=NaN;
    R.Ev(k).ttestEvent(1:Nneurons,1)=NaN;
    R.Ev(k).RespDirPre(1:Nneurons,1)=NaN;
    R.Ev(k).RespDirPost(1:Nneurons,1)=NaN;
    R.Ev(k).RespDirEvent(1:Nneurons,1)=NaN;
    R.Ev(k).ActDirPre(1:Nneurons,1)=NaN;
    R.Ev(k).ActDirPost(1:Nneurons,1)=NaN;
    R.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
    R.Ev(k).MaxVal(1:Nneurons,1)=NaN;
    R.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 R.Structure(NN)~=0 %to avoid analyzing excluded neurons

            ev2=NaN(size(RAW(i).Eint{1,3},1),2);
            [PSR1iti,N1iti]=MakePSR0int(RAW(i).Nrast(j),RAW(i).Eint{1,1},RAW(i).Eint{1,2},{2},{1});% makes trial by trial rasters for baseline based on ITI
            [PSR1trial,N1trial]=MakePSR0int(RAW(i).Nrast(j),RAW(i).Eint{1,3},RAW(i).Eint{1,4},{2},{2});% makes trial by trial rasters for each trial based on time interval

            for l=1:length(RAW(i).Eint{1,3}) %loops through trials
                ev2(l,1)=sum(PSR1iti{l}<RAW(i).Eint{1,2}(l) & PSR1iti{l}>RAW(i).Eint{1,1}(l))/(RAW(i).Eint{1,2}(l)-RAW(i).Eint{1,1}(l)-20);
                ev2(l,2)=sum(PSR1trial{l}<RAW(i).Eint{1,4}(l) & PSR1trial{l}>RAW(i).Eint{1,3}(l))/(RAW(i).Eint{1,4}(l)-RAW(i).Eint{1,3}(l));
            end
             
             R.WF(NN,:)=RAW(i).waveforms(j,:);   
             R.WFanalysis(NN,1)=mean(ev2(:,1)); % mean FR during Baseline in the middle of each ITI (first and last 10s excluded)
             proportionISI=0;
             for k=3:1:size(RAW(i).Nrast{j,1},1) % loop thru timestamps waveform
                 ISIn=RAW(i).Nrast{j,1}(k)-RAW(i).Nrast{j,1}(k-1); 
                 ISIo=RAW(i).Nrast{j,1}(k-1)-RAW(i).Nrast{j,1}(k-2);
                 calcul(k-2)=2*abs(ISIn-ISIo)/(ISIn+ISIo);
                 if ISIo>2
                     proportionISI=proportionISI+1;        
                 end
             end
             R.CV2(NN,1)=mean(calcul);
             R.CV2(NN,2)=100*proportionISI/(size(RAW(i).Nrast{j,1},1)-1);% of ISI>2s
         

            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
                
                R.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
                Tbase=prewin(k,1):BSIZE:prewin(k,2);
                    [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)
%                     if strcmp(Erefnames(k),'BaselineITI') 
%                         [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},DuraITI,{2});% makes trial by trail rasters for BL based on ITI. PSR1 is a cell(neurons, trials)
%                     else
                    [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
                    %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
                        
                        % calculate MeanZ from PSTH before smoothing
                            if sum(PTH0,2)~=0
                                unsmoothPSTHz(1,1:length(Tm))=normalize(PTH1,PTH0,0);
                                R.Ev(k).rawMeanz(NN,1)=nanmean(unsmoothPSTHz(1,Tm>EventWin(k,1) & Tm<EventWin(k,2)),2);
                                R.Ev(k).rawMeanzPre(NN,1)=nanmean(unsmoothPSTHz(1,Tm>PreEventWin(k,1) & Tm<PreEventWin(k,2)),2);
                            else
                                R.Ev(k).rawMeanz(NN,1)=NaN;
                                R.Ev(k).rawMeanzPre(NN,1)=NaN;
                            end


                        PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
                        PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';
                        
                        %------------- Fills the R.Ev(k) fields --------------
                        %R.Ev(k).BW(NN,1)=BW1;
                        R.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
                        R.Ev(k).PSTHrawBL(NN,1:length(Tbase))=PTH0;
                        R.Ev(k).Meanraw(NN,1)=nanmean(R.Ev(k).PSTHraw(NN,Tm>PostEventWin(k,1) & Tm<PostEventWin(k,2)),2);
                        if sum(PTH0,2)~=0
                            R.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
                            R.Ev(k).Meanz(NN,1)=nanmean(R.Ev(k).PSTHz(NN,Tm>EventWin(k,1) & Tm<EventWin(k,2)),2);
                            R.Ev(k).MeanzPRE(NN,1)=nanmean(R.Ev(k).PSTHz(NN,Tm>PreEventWin(k,1) & Tm<PreEventWin(k,2)),2);
                        else
                            R.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
                            R.Ev(k).Meanz(NN,1)=NaN;
                            R.Ev(k).MeanzPRE(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).PreEvent=NaN(size(RAW(i).Erast{EvInd},1),1);
                        ev(k).PostEvent=NaN(size(RAW(i).Erast{EvInd},1),1);
                        ev(k).Event=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).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));
                            ev(k).Event(m)=sum(PSR2{m}<EventWin(k,2) & PSR2{m}>EventWin(k,1))/(EventWin(k,2)-EventWin(k,1));
                        end
                              
                        %---------------------------Response detection w/ SignRank on pre/post windows---------------------------

                    if  ~isempty(EvInd) && R.Ev(k).NumberTrials(NN,1)>= MinNumTrials && nanmean(ev(k).pre(:,1),1)>0.5 %avoid analyzing sessions that do not have enough trials / exclude low firing neurons
                        
                        [~,R.Ev(k).ttestPreEvent(NN,1)]=ttest(ev(k).pre, ev(k).PreEvent); %Signrank used here because it is a dependant sample test
                        if R.Ev(k).ttestPreEvent(NN,1)<PStat
                            R.Ev(k).RespDirPre(NN,1)=sign(mean(ev(k).PreEvent)-mean(ev(k).pre));                           
                        else R.Ev(k).RespDirPre(NN,1)=0;
                        end
                        
                        [~,R.Ev(k).ttestPostEvent(NN,1)]=ttest(ev(k).pre, ev(k).PostEvent); %Signrank used here because it is a dependant sample test
                        if R.Ev(k).ttestPostEvent(NN,1)<PStat
                            R.Ev(k).RespDirPost(NN,1)=sign(mean(ev(k).PostEvent)-mean(ev(k).pre));
                        else R.Ev(k).RespDirPost(NN,1)=0;
                        end
                        
                        [~,R.Ev(k).ttestEvent(NN,1)]=ttest(ev(k).pre, ev(k).Event); %Signrank used here because it is a dependant sample test
                        if R.Ev(k).ttestEvent(NN,1)<PStat
                            R.Ev(k).RespDirEvent(NN,1)=sign(mean(ev(k).Event)-mean(ev(k).pre));
                        else R.Ev(k).RespDirEvent(NN,1)=0;
                        end
                            
                        if R.Ev(k).RespDirPre(NN,1)>0 || R.Ev(k).RespDirPost(NN,1) >0
                            Search=find(Tm>=PreEventWin(k,1) & Tm<=PostEventWin(k,2));
                            [R.Ev(k).MaxVal(NN,1),MaxInd]=max(R.Ev(k).PSTHraw(NN,Search));
                            R.Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                        elseif R.Ev(k).RespDirPre(NN,1)<0 || R.Ev(k).RespDirPost(NN,1) <0
                            Search=find(Tm>=PreEventWin(k,1) & Tm<=PostEventWin(k,2));
                            [R.Ev(k).MaxVal(NN,1),MaxInd]=min(R.Ev(k).PSTHraw(NN,Search));
                            R.Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                        end
                        
                        TRNeventPOST(NN,k)=R.Ev(k).RespDirPost(NN);
                        TRNeventPRE(NN,k)=R.Ev(k).RespDirPre(NN);
                        
                        %---------------------------Response detection w/ RESDETECT08 on running windows---------------------------

                    end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
                end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN OR low firing during baseline
            end %Events
           
            %%----------------------- STAT -----------------------
  
            fprintf('Neuron ID # %d\n',NN);
        elseif R.Structure(NN)~=0
        end %exclusion: IF R.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)

TRNevent=[TRNeventPRE(:,1:6) TRNeventPOST(:,1:5) TRNeventPOST(:,7)]; 
R.TRN(:,1)=sum(abs(TRNevent(NN,:)),2,'omitnan'); % if >0 : Task related responsive neurons (include excitation and inhibition)
R.TRN(:,2)=sum(TRNevent(NN,:),2,'omitnan'); % if >0 : TRN excited. if < 0: TRN inhibited. may loose ambivalent neurons (both excited and inhib)
         
if SAVE_FLAG
 save('Rextendedtraining.mat','R')
end

toc