3Tasks_paper/Scripts/D_Figure4.m

%% Plotting Heatmaps and average PSTH for sequence of event in DT5, FR5 and FS5
% edited on june 2019 to only include sucrose rats in DT5 data set
% do average first/last 3 session on PSTH instead of all sessions

clc;clear all; close all;


task_name={'DT5','FR5','FS5'};
task_Dataset={'R_DT5','R_FR5','R_FS5'};
BSIZE=0.01; %Do not change
Dura=[-25 20]; Tm=Dura(1):BSIZE:Dura(2);

Eventnames{1}={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
Eventnames{2}={'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
Eventnames{3}={'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};

%% --- --- Plotting color-coded and average PSTHs --- ---
BrainRegion=[10 20]; %10 for DLS, 20 for DMS

% ----------- COLORS ---------
myblue=[0 113/255 188/255];
Clim=[-2 2];%ColorMap
Xevent=[-0.25 0.25];
Yaxis=[-1 3];Xaxis=[1 357];
Ishow=find(Tm>=Xevent(1) & Tm<=Xevent(2));

% for stat
Ishow1(1,:)=find(Tm>=Xevent(1) & Tm<=0);
Ishow1(2,:)=find(Tm>=0 & Tm<=Xevent(2));
Combined_Table_allTask=[];allTask_ID=[];

%% ------------------------------------Main program ------------------------------------------------------------------------------------------------------------------------
figure;
for task=1:length(task_Dataset)
    %generate concat file
    clear R DLSselectionTRN DMSselectionTRN Celltype SesCelltype
    load(task_Dataset{task})
    DLSconcat_ACQ_first3=[];DMSconcat_ACQ_first3=[];DLSconcat_ACQ_last3=[]; DMSconcat_ACQ_last3=[];DLSconcat_ACQ_all=[];DMSconcat_ACQ_all=[];
    for k=1:size(R.Ses,2)
        DLSselectionTRN=R.Ses(k).Coord(:,4)==BrainRegion(1) & R.Ses(k).TRN(:,1)~=0  & SesCelltype(k).Celltype(:,1)==1 ;
        DMSselectionTRN=R.Ses(k).Coord(:,4)==BrainRegion(2) & R.Ses(k).TRN(:,1)~=0  & SesCelltype(k).Celltype(:,1)==1;
        NbN{task}(k,1)=sum(DLSselectionTRN)*100/sum(R.Ses(k).Coord(:,4)==BrainRegion(1) & SesCelltype(k).Celltype(:,1)==1);
        NbN{task}(k,2)=100-NbN{1}(k,1);
        NbN{task}(k,3)=sum(DMSselectionTRN)*100/sum(R.Ses(k).Coord(:,4)==BrainRegion(2) & SesCelltype(k).Celltype(:,1)==1);
        NbN{task}(k,4)=100-NbN{1}(k,3);
        
        DLSconcat_ACQ=[];DMSconcat_ACQ=[];
        if task==1
            DLSconcat_ACQ = cat(2,DLSconcat_ACQ,R.Ses(k).Ev(7).PSTHz(DLSselectionTRN,Ishow));
            DMSconcat_ACQ = cat(2,DMSconcat_ACQ,R.Ses(k).Ev(7).PSTHz(DMSselectionTRN,Ishow));
        end
        for i=1:6
            DLSconcat_ACQ = cat(2,DLSconcat_ACQ,R.Ses(k).Ev(i).PSTHz(DLSselectionTRN,Ishow));
            DMSconcat_ACQ = cat(2,DMSconcat_ACQ,R.Ses(k).Ev(i).PSTHz(DMSselectionTRN,Ishow));
        end
        
        if k<4 %first 3 DT5 sessions
            DLSconcat_ACQ_first3=cat(1,DLSconcat_ACQ_first3,DLSconcat_ACQ);
            DMSconcat_ACQ_first3=cat(1,DMSconcat_ACQ_first3,DMSconcat_ACQ);
        end
        if k>(length(R.Ses)-3)
            DLSconcat_ACQ_last3=cat(1,DLSconcat_ACQ_last3,DLSconcat_ACQ);
            DMSconcat_ACQ_last3=cat(1,DMSconcat_ACQ_last3,DMSconcat_ACQ);
        end
        DLSconcat_ACQ_all=cat(1,DLSconcat_ACQ_all,DLSconcat_ACQ);
        DMSconcat_ACQ_all=cat(1,DMSconcat_ACQ_all,DMSconcat_ACQ);
        time=[1:1:size(DLSconcat_ACQ,2)];
        
        MeanDLS_ACQ{task}(k,:)=nanmean(DLSconcat_ACQ,1);
        MeanDMS_ACQ{task}(k,:)=nanmean(DMSconcat_ACQ,1);
        DLSsem1{task}(k,:)=nansem(DLSconcat_ACQ(:,time),1);
        DMSsem1{task}(k,:)=nansem(DMSconcat_ACQ(:,time),1);
        DLSnn{task}(k,1)=size(DLSconcat_ACQ,1);
        DMSnn{task}(k,1)=size(DMSconcat_ACQ,1);
    end
    MeanDLS_ACQ_first3=nanmean(DLSconcat_ACQ_first3,1);
    MeanDMS_ACQ_first3=nanmean(DMSconcat_ACQ_first3,1);
    DLSsem1_first3=nansem(DLSconcat_ACQ_first3,1);
    DMSsem1_first3=nansem(DMSconcat_ACQ_first3,1);
    
    MeanDLS_ACQ_last3=nanmean(DLSconcat_ACQ_last3,1);
    MeanDMS_ACQ_last3=nanmean(DMSconcat_ACQ_last3,1);
    DLSsem1_last3=nansem(DLSconcat_ACQ_last3,1);
    DMSsem1_last3=nansem(DMSconcat_ACQ_last3,1);
    
    MeanDLS_ACQ_all=nanmean(DLSconcat_ACQ_all,1);
    MeanDMS_ACQ_all=nanmean(DMSconcat_ACQ_all,1);
    DLSsem1_all=nansem(DLSconcat_ACQ_all,1);
    DMSsem1_all=nansem(DMSconcat_ACQ_all,1);
    
    %% table for stat
    Combined_Session=[]; Combined_Table=[];Combined_region=[];
    for k=1:size(R.Ses,2)
        DSconcat_ACQ=[];Table=[];
        for l=1:6
            DSconcat_ACQ = cat(2,DSconcat_ACQ,R.Ses(k).Ev(l).PSTHz(:,Ishow));
        end
        selection=sum(isnan(DSconcat_ACQ),2)==0 & R.Ses(k).Coord(:,4)~=0 & R.Ses(k).TRN(:,1)~=0 & SesCelltype(k).Celltype(:,1)==1;
        for l=1:6
            for m=1:2
                Table=cat(2,Table,nanmean(R.Ses(k).Ev(l).PSTHz(selection,Ishow1(m,:)),2));
            end
        end
        Combined_Table=cat(1,Combined_Table,Table);
        Combined_Session=cat(1,Combined_Session,k*ones(size(Table,1),1));
        Combined_region=cat(1,Combined_region,R.Ses(k).Coord(selection,4));
    end
    Combined_Table=[Combined_Session Combined_region Combined_Table];
    Combined_Table_allTask=cat(1,Combined_Table_allTask,Combined_Table);
    allTask_ID=cat(1,allTask_ID, task*ones(length(Combined_Table),1));
    TableActivity=cat(2,allTask_ID,Combined_Table_allTask);
    %% Heatmaps
    if task>1
        Xaxis=[1 306];
    end
    
    figure(1)
    subplot(3,3,1+3*(task-1))
    imagesc(time,[1 size(MeanDLS_ACQ{task},1)],MeanDLS_ACQ{task},Clim); %colorbar;axis tight;
    colormap(parula);
    hold on
    plot([51 51],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    plot([102 102],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    plot([153 153],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    plot([204 204],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    plot([255 255],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    plot([306 306],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    hold off
    set(gca,'XTick',0:25.5:size(MeanDLS_ACQ{task},2));
    set(gca,'xticklabel',Eventnames{task})
    ylabel('sessions');
    
    
    subplot(3,3,2+3*(task-1))
    imagesc(time,[1 size(MeanDMS_ACQ{task},1)],MeanDMS_ACQ{task},Clim); %colorbar;axis tight;
    colormap(parula);
    hold on
    plot([51 51],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    plot([102 102],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    plot([153 153],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    plot([204 204],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    plot([255 255],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    plot([306 306],[size(MeanDLS_ACQ{task},1)+0.5 0.5],'k')
    set(gca,'XTick',0:25.5:size(MeanDLS_ACQ{task},2));
    set(gca,'xticklabel',Eventnames{task})
    hold off
    
    subplot(3,3,3+3*(task-1))
    DLSpsth=MeanDLS_ACQ_all;
    DLSsem=DLSsem1_all;
    DLSup=DLSpsth+DLSsem;
    DLSdown=DLSpsth-DLSsem;
    hold on;
    patch([time,time(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
    g1=plot(time,DLSpsth,'Color',myblue,'linewidth',1.5);
    
    DMSpsth=MeanDMS_ACQ_all;
    DMSsem=DMSsem1_all;
    DMSup=DMSpsth+DMSsem;
    DMSdown=DMSpsth-DMSsem;
    hold on;
    patch([time,time(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
    g2=plot(time,DMSpsth,'r','linewidth',1.5);
    plot([0 357], [0 0],'LineStyle',':','Color','k');
    plot([51 51],[-1 3],'LineStyle',':','Color','k')
    plot([102 102],[-1 3],'LineStyle',':','Color','k')
    plot([153 153],[-1 3],'LineStyle',':','Color','k')
    plot([204 204],[-1 3],'LineStyle',':','Color','k')
    plot([255 255],[-1 3],'LineStyle',':','Color','k')
    plot([306 306],[-1 3],'LineStyle',':','Color','k')
    axis([Xaxis,Yaxis]);
    set(gca,'XTick',0:25.5:size(MeanDLS_ACQ,2));
    set(gca,'xticklabel',Eventnames{task})
    ylabel('z-score');
    
    trace=cat(2,g1,g2);
    label{1}=['DLS'];label{2}=['DMS'];
    legend(trace,label,'location','northeast');
    hold off
    
    %% average PSTH
    figure(2)
    subplot(6,1,1+2*(task-1))
    DLSpsth=MeanDLS_ACQ_first3;
    DLSsem=DLSsem1_first3;
    DLSup=DLSpsth+DLSsem;
    DLSdown=DLSpsth-DLSsem;
    hold on;
    patch([time,time(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
    g1=plot(time,DLSpsth,'Color',myblue,'linewidth',1.5);
    
    DMSpsth=MeanDMS_ACQ_first3;
    DMSsem=DMSsem1_first3;
    DMSup=DMSpsth+DMSsem;
    DMSdown=DMSpsth-DMSsem;
    hold on;
    patch([time,time(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
    g2=plot(time,DMSpsth,'r','linewidth',1.5);
    plot([0 357], [0 0],'LineStyle',':','Color','k');
    plot([51 51],[-1 3],'LineStyle',':','Color','k')
    plot([102 102],[-1 3],'LineStyle',':','Color','k')
    plot([153 153],[-1 3],'LineStyle',':','Color','k')
    plot([204 204],[-1 3],'LineStyle',':','Color','k')
    plot([255 255],[-1 3],'LineStyle',':','Color','k')
    plot([306 306],[-1 3],'LineStyle',':','Color','k')
    axis([Xaxis,Yaxis]);
    set(gca,'XTick',0:25.5:size(MeanDLS_ACQ,2));
    set(gca, 'XTickLabel', [])
    ylabel('z-score');
    
    trace=cat(2,g1,g2);
    label{1}=['DLS'];label{2}=['DMS'];
    legend(trace,label,'location','northeast');
    hold off
    
    subplot(6,1,2+2*(task-1))
    DLSpsth=MeanDLS_ACQ_last3;
    DLSsem=DLSsem1_last3;
    DLSup=DLSpsth+DLSsem;
    DLSdown=DLSpsth-DLSsem;
    hold on;
    patch([time,time(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
    g1=plot(time,DLSpsth,'Color',myblue,'linewidth',1.5);
    
    DMSpsth=MeanDMS_ACQ_last3;
    DMSsem=DMSsem1_last3;
    DMSup=DMSpsth+DMSsem;
    DMSdown=DMSpsth-DMSsem;
    hold on;
    patch([time,time(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
    g2=plot(time,DMSpsth,'r','linewidth',1.5);
    plot([0 357], [0 0],'LineStyle',':','Color','k');
    plot([51 51],[-1 3],'LineStyle',':','Color','k')
    plot([102 102],[-1 3],'LineStyle',':','Color','k')
    plot([153 153],[-1 3],'LineStyle',':','Color','k')
    plot([204 204],[-1 3],'LineStyle',':','Color','k')
    plot([255 255],[-1 3],'LineStyle',':','Color','k')
    plot([306 306],[-1 3],'LineStyle',':','Color','k')
    axis([Xaxis,Yaxis]);
    set(gca,'XTick',0:25.5:size(MeanDLS_ACQ,2));
    set(gca,'xticklabel',Eventnames{task})
    ylabel('z-score');
    
end

%% statistiques
nbevent=[1:12];

stat(1).condition{1}='TRN DT5';
stat(2).condition{1}='TRN FR5';
stat(3).condition{1}='TRN FS5';
filename = 'ResultStat_Activity_Fig4.xlsx'; 
for i=1:3
    selection=TableActivity(:,1)==i;
    tableSTAT=array2table(TableActivity(selection,2:end),'VariableNames',{'session','region','win1','win2','win3','win4','win5','win6','win7','win8','win9','win10','win11','win12'});
    tableSTAT.session=categorical(tableSTAT.session);
    tableSTAT.region=categorical(tableSTAT.region);
    
    rm_reg = fitrm(tableSTAT,'win1-win12~region*session','WithinDesign',nbevent');
    stat(i).ranovatbl_event = ranova(rm_reg);
    stat(i).Btw_ranovatbl_event = anova(rm_reg);

    writetable(stat(i).ranovatbl_event,filename,'Sheet',i,'Range','A1')
    writetable(stat(i).Btw_ranovatbl_event,filename,'Sheet',i,'Range','J1')
end

stat(4).condition={'3 tasks'};
selection=TableActivity(:,2)<8; % Comparison tasks, sessions 1-7
tableSTAT=array2table(TableActivity(selection,:),'VariableNames',{'task','session','region','win1','win2','win3','win4','win5','win6','win7','win8','win9','win10','win11','win12'});
tableSTAT.session=categorical(tableSTAT.session);
tableSTAT.region=categorical(tableSTAT.region);
tableSTAT.task=categorical(tableSTAT.task);

rm_reg = fitrm(tableSTAT,'win1-win12~task*session*region','WithinDesign',nbevent');
stat(4).ranovatbl_event = ranova(rm_reg);
stat(4).Btw_ranovatbl_event = anova(rm_reg);
stat(4).posthoc = multcompare(rm_reg,'task', 'by', 'region');

writetable(stat(4).ranovatbl_event,filename,'Sheet',4,'Range','A1')
writetable(stat(4).Btw_ranovatbl_event,filename,'Sheet',4,'Range','J1')
writetable(stat(4).posthoc,filename,'Sheet',4,'Range','R1')