auksztulewicz_et_al_2023_bmcbiol/scripts/Figure3A_MUA.m

clear all

new_fs=150;
t_ax=[-.1:1/new_fs:.4];
baseline=[-.05 0.0133];
timewin=[-.05 .25];
t_ax_plot=t_ax(min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))));
badchan=48;
t_ax_short=t_ax_plot;

load('C:\Users\Ryzard Auksztulewicz\Desktop\omissions_bmcbiol_repo\data\ryszard_preprocessed_20200413_AMUA.mat')

for p=[5 8:12]

    rates_exp=[];
    mean_response=[];
    seg_length=[];
    cols=cell(0);
    linestyles=cell(0);
    legendlabels=cell(0);

    for j=1:3

        data_present=length(preceding_all_amua{p,j});
        

        if data_present>0


            notrls=size(preceding_all_amua{p,j},3);

            for k=1:size(preceding_all_amua{p,j},3)
                for l=1:size(preceding_all_amua{p,j},1)
                    for_detrend=[mean(preceding_all_amua{p,j}(l,min(find(t_ax>=baseline(1))):max(find(t_ax<=baseline(2))),k)) mean(preceding_all_amua{p,j}(l,min(find(t_ax>=timewin(2)+baseline(1))):max(find(t_ax<=timewin(2))),k))];
                    for_detrend=for_detrend-mean(for_detrend);
                    preceding_all_amua_dt{p,j}(l,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),k)=preceding_all_amua{p,j}(l,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),k)-linspace(for_detrend(1),for_detrend(2),length(find(t_ax>=timewin(1)&t_ax<=timewin(2))));
                end
            end
            for k=1:size(preceding_all_amua{p,j},3)
                for l=1:size(preceding_all_amua{p,j},1)
                    for_detrend=[mean(omissions_all_amua{p,j}(l,min(find(t_ax>=baseline(1))):max(find(t_ax<=baseline(2))),k)) mean(omissions_all_amua{p,j}(l,min(find(t_ax>=timewin(2)+baseline(1))):max(find(t_ax<=timewin(2))),k))];
                    for_detrend=for_detrend-mean(for_detrend);
                    omissions_all_amua_dt{p,j}(l,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),k)=omissions_all_amua{p,j}(l,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),k)-linspace(for_detrend(1),for_detrend(2),length(find(t_ax>=timewin(1)&t_ax<=timewin(2))));
                end
            end

            preceding_all_amua_bc{p,j}=preceding_all_amua_dt{p,j}-repmat(nanmean(preceding_all_amua_dt{p,j}(:,min(find(t_ax>=baseline(1))):max(find(t_ax<=baseline(2))),:),2),[1 size(preceding_all_amua_dt{p,j},2) 1]);
            omissions_all_amua_bc{p,j}=omissions_all_amua_dt{p,j}-repmat(nanmean(omissions_all_amua_dt{p,j}(:,min(find(t_ax>=baseline(1))):max(find(t_ax<=baseline(2))),:),2),[1 size(omissions_all_amua_dt{p,j},2) 1]);

            preceding_all_amua_bc{p,j}=preceding_all_amua_bc{p,j}./repmat(nanstd(preceding_all_amua_bc{p,j}(:,min(find(t_ax>=baseline(1))):max(find(t_ax<=baseline(2))),:),[],2),[1 size(preceding_all_amua_bc{p,j},2) 1]);
            omissions_all_amua_bc{p,j}=omissions_all_amua_bc{p,j}./repmat(nanstd(omissions_all_amua_bc{p,j}(:,min(find(t_ax>=baseline(1))):max(find(t_ax<=baseline(2))),:),[],2),[1 size(omissions_all_amua_bc{p,j},2) 1]);
        end
    end

    mean_evoked=[];
    for j=1:3
        if length(preceding_all_amua{p,j})>0
            mean_evoked=[mean_evoked nanmean(preceding_all_amua_bc{p,j}(:,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),:),3)];
        end
    end
    totalrange=range(mean_evoked,2);
    for j=1:3
        if length(preceding_all_amua{p,j})>0
            preceding_all_amua_bc{p,j}=preceding_all_amua_bc{p,j}./repmat(totalrange,[1 size(preceding_all_amua_bc{p,j},2) size(preceding_all_amua_bc{p,j},3)]);
            omissions_all_amua_bc{p,j}=omissions_all_amua_bc{p,j}./repmat(totalrange,[1 size(omissions_all_amua_bc{p,j},2) size(omissions_all_amua_bc{p,j},3)]);
        end
    end

    allstd=[];
    for j=1:size(preceding_all_amua_bc,2)
        try
            allstd=[allstd; squeeze(mean(std(omissions_all_amua_bc{p,j}(setdiff(1:64,badchan),find(t_ax>=timewin(1)&t_ax<=timewin(2)),:),[],2),1))];
        catch
        end
    end
    stdcutoff=median(allstd)+3*std(allstd);
    for j=1:size(preceding_all_amua_bc,2)
        try
            badtrls=find(squeeze(mean(std(omissions_all_amua_bc{p,j}(setdiff(1:64,badchan),find(t_ax>=timewin(1)&t_ax<=timewin(2)),:),[],2),1))>stdcutoff);
            omissions_all_amua_bc{p,j}(:,:,badtrls)=NaN;
        catch
        end
    end

    allstd=[];
    for j=1:size(preceding_all_amua_bc,2)
        try
            allstd=[allstd; squeeze(mean(std(preceding_all_amua_bc{p,j}(setdiff(1:64,badchan),find(t_ax>=timewin(1)&t_ax<=timewin(2)),:),[],2),1))];
        catch
        end
    end
    stdcutoff=median(allstd)+3*std(allstd);
    for j=1:size(preceding_all_amua_bc,2)
        try
            badtrls=find(squeeze(mean(std(preceding_all_amua_bc{p,j}(setdiff(1:64,badchan),find(t_ax>=timewin(1)&t_ax<=timewin(2)),:),[],2),1))>stdcutoff);
            preceding_all_amua_bc{p,j}(:,:,badtrls)=NaN;
        catch
        end
    end

    for j=1:3
        if length(preceding_all_amua{p,j})>0


            mean_response=[mean_response nanmean(preceding_all_amua_bc{p,j}(:,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),:),3) nanmean(omissions_all_amua_bc{p,j}(:,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),:),3)];
            rates_exp=[rates_exp j+1];
            seg_length=[seg_length size(nanmean(preceding_all_amua_bc{p,j}(:,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),:),3),2) size(nanmean(omissions_all_amua_bc{p,j}(:,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),:),3),2)];
        end
    end

    mean_response(badchan,:)=0;




    for c=1:64
        alltrls=[];
        for k=1:size(preceding_all_amua_bc,2)
            try
                alltrls=[alltrls nanmean(squeeze(preceding_all_amua_bc{p,k}(c,find(t_ax>=0&t_ax<=timewin(2)),:)),1)];
            catch
            end
        end
        [h,p_prec(c),ci,stats]=ttest(alltrls);
    end
    for c=1:64
        alltrls=[];
        for k=1:size(omissions_all_amua_bc,2)
            try
                alltrls=[alltrls nanmean(squeeze(omissions_all_amua_bc{p,k}(c,find(t_ax>=0&t_ax<=timewin(2)),:)),1)];
            catch
            end
        end
        [h,p_omi(c),ci,stats]=ttest(alltrls);
    end

    p_sort=sort([p_prec p_omi]);
    for i=1:length(p_sort)
        if p_sort(i)<=.05*i/length(p_sort)
            p_fdr(i)=1;
        else
            p_fdr(i)=0;
        end
    end
    p_fdr_thresh=p_sort(max(find(p_fdr==1)));
    if length(p_fdr_thresh)==0
        p_fdr_thresh=0;
    end
    selchans=setdiff(1:64,badchan);

    try
        amps_all=[];
        lats_all=[];
        for c=selchans
            mean_conds=reshape(mean_response(c,:),[45 size(mean_response,2)/45]);
            if mean(sign(mean(mean_conds)))<0
                [p c]
                mean_conds=-mean_conds;
                mean_response(c,:)=-mean_response(c,:);
            end
             for k=1:size(mean_response,2)/45
                [amps,lats]=findpeaks(mean_conds(:,k));
                lats_all(c,k)=t_ax_short(lats(find(abs(amps)==max(abs(amps)))));
                amps_all(c,k)=amps(find(abs(amps)==max(abs(amps))));
             end
        end
        lats_all=lats_all(selchans,:);
        amps_all=amps_all(selchans,:);

        mean_conds_all{p}=mean_response;
        mean_lats_all{p}=lats_all;
        mean_amps_all{p}=amps_all;
        selchans_all{p}=selchans;
        rates_all{p}=rates_exp;
    catch
        mean_conds_all{p}=0;
        mean_lats_all{p}=0;
        mean_amps_all{p}=0;
        selchans_all{p}=0;
        rates_all{p}=0;
        p
    end

end


choose_rates=[2 3 4];
anova_lats=[];
anova_amps=[];
anova_rats=[];
anova_rate=[];
anova_stim=[];
anova_chan=[];
ppoi=[5 8:12];

colscols=[0 0 .5; 0 0 1; .5 .5 1];

for i=ppoi
    if length(rates_all{i})>0 && length([find(rates_all{i}==choose_rates(1)) find(rates_all{i}==choose_rates(2))])==2

        selchoi=find(median(mean_amps_all{i}(:,[2 4 6]),2)<0);
        mean_lats_all{i}(selchoi,:)=NaN;
        mean_amps_all{i}(selchoi,:)=NaN;
        for r=1:length(choose_rates)
            rateind=find(rates_all{i}==choose_rates(r));
            nchans=length(find(~isnan(mean_lats_all{i}(:,(rateind-1)*2+1))));
            anova_lats=[anova_lats; mean_lats_all{i}(find(~isnan(mean_lats_all{i}(:,(rateind-1)*2+1))),(rateind-1)*2+1)];
            anova_amps=[anova_amps; mean_amps_all{i}(find(~isnan(mean_amps_all{i}(:,(rateind-1)*2+1))),(rateind-1)*2+1)];
            anova_rats=[anova_rats; ones(nchans,1)*i];
            anova_rate=[anova_rate; ones(nchans,1)*r];
            anova_stim=[anova_stim; ones(nchans,1)*1];
            anova_chan=[anova_chan; (1:nchans)'];
            anova_lats=[anova_lats; mean_lats_all{i}(find(~isnan(mean_lats_all{i}(:,(rateind-1)*2+2))),(rateind-1)*2+2)];
            anova_amps=[anova_amps; mean_amps_all{i}(find(~isnan(mean_amps_all{i}(:,(rateind-1)*2+2))),(rateind-1)*2+2)];
            anova_rats=[anova_rats; ones(nchans,1)*i];
            anova_rate=[anova_rate; ones(nchans,1)*r];
            anova_stim=[anova_stim; ones(nchans,1)*2];
            anova_chan=[anova_chan; (1:nchans)'];
        end
    end
end
M = zeros(4,4);
M(3,4) = 1;


anovan(anova_amps(~isnan(anova_amps)),{anova_stim(~isnan(anova_amps)) anova_rate(~isnan(anova_amps))  anova_rats(~isnan(anova_amps))},'varnames',{'stim' 'rate'  'rat'},'random',[3],'model','full','display','on')
anovan(anova_lats(~isnan(anova_amps)),{anova_stim(~isnan(anova_amps)) anova_rate(~isnan(anova_amps))  anova_rats(~isnan(anova_amps))},'varnames',{'stim' 'rate'  'rat'},'random',[3],'model','full','display','on')

figure;
subplot(1,4,2)
hold on
unique_rats=unique(anova_rats);
for i=1:length(unique_rats)
    for j=1:length(choose_rates)
        scatter(mean(anova_lats(find(anova_rats==unique_rats(i)&anova_stim==1&anova_rate==choose_rates(j)-1))),mean(anova_amps(find(anova_rats==unique_rats(i)&anova_stim==1&anova_rate==choose_rates(j)-1))),25,colscols(j,:),'filled')
        scatter(mean(anova_lats(find(anova_rats==unique_rats(i)&anova_stim==2&anova_rate==choose_rates(j)-1))),mean(anova_amps(find(anova_rats==unique_rats(i)&anova_stim==2&anova_rate==choose_rates(j)-1))),25,colscols(j,:))
        line([mean(anova_lats(find(anova_rats==unique_rats(i)&anova_stim==1&anova_rate==choose_rates(j)-1))) mean(anova_lats(find(anova_rats==unique_rats(i)&anova_stim==2&anova_rate==choose_rates(j)-1)))],[mean(anova_amps(find(anova_rats==unique_rats(i)&anova_stim==1&anova_rate==choose_rates(j)-1))) mean(anova_amps(find(anova_rats==unique_rats(i)&anova_stim==2&anova_rate==choose_rates(j)-1)))],'color',colscols(j,:))
    end
end
xlim([timewin(1) timewin(2)])
xlabel('peak latency (s)')
ylabel('peak amplitude')
line([-.2 1],[0 0],'linestyle','--','color',[.5 .5 .5])
line([0 0],[-.2 1],'linestyle','--','color',[.5 .5 .5])

subplot(1,4,3)
hold on
for i=1:length(choose_rates)
    plotmean=mean(anova_lats(find(anova_stim==2&anova_rate==choose_rates(i)-1)));
    plotsem=std(anova_lats(find(anova_stim==2&anova_rate==choose_rates(i)-1)))/sqrt(length(anova_lats(find(anova_stim==2&anova_rate==choose_rates(i)-1))));
    bar(3*choose_rates(i)+.5,plotmean,'facecolor','white','edgecolor',colscols(i,:))
    line([3*choose_rates(i)+.5 3*choose_rates(i)+.5],[plotmean-plotsem plotmean+plotsem],'color',colscols(i,:))
    plotmean=mean(anova_lats(find(anova_stim==1&anova_rate==choose_rates(i)-1)));
    plotsem=std(anova_lats(find(anova_stim==1&anova_rate==choose_rates(i)-1)))/sqrt(length(anova_lats(find(anova_stim==2&anova_rate==choose_rates(i)-1))));
    bar(3*choose_rates(i)-.5,plotmean,'facecolor',colscols(i,:))
    line([3*choose_rates(i)-.5 3*choose_rates(i)-.5],[plotmean-plotsem plotmean+plotsem],'color',colscols(i,:))
end
ylabel('peak latency (s)')
xlim([4 14])
xticks(6:3:12)
xticklabels(2:4)
xlabel('rate (Hz)')
ylim([0 .16])

subplot(1,4,4)
hold on
for i=1:length(choose_rates)
    plotmean=mean(anova_amps(find(anova_stim==2&anova_rate==choose_rates(i)-1)));
    plotsem=std(anova_amps(find(anova_stim==2&anova_rate==choose_rates(i)-1)))/sqrt(length(anova_lats(find(anova_stim==2&anova_rate==choose_rates(i)-1))));
    bar(3*choose_rates(i)+.5,plotmean,'facecolor','white','edgecolor',colscols(i,:))
    line([3*choose_rates(i)+.5 3*choose_rates(i)+.5],[plotmean-plotsem plotmean+plotsem],'color',colscols(i,:))
    plotmean=mean(anova_amps(find(anova_stim==1&anova_rate==choose_rates(i)-1)));
    plotsem=std(anova_amps(find(anova_stim==1&anova_rate==choose_rates(i)-1)))/sqrt(length(anova_lats(find(anova_stim==2&anova_rate==choose_rates(i)-1))));
    bar(3*choose_rates(i)-.5,plotmean,'facecolor',colscols(i,:))
    line([3*choose_rates(i)-.5 3*choose_rates(i)-.5],[plotmean-plotsem plotmean+plotsem],'color',colscols(i,:))
end
ylabel('peak amplitude')
xlim([4 14])
xticks(6:3:12)
xticklabels(2:4)
xlabel('rate (Hz)')
ylim([0 .8])

clear all_ps
subplot(1,4,1)
plot_data=[mean_conds_all{5}(selchans_all{5},:); mean_conds_all{8}(selchans_all{8},:); mean_conds_all{9}(selchans_all{9},:); mean_conds_all{10}(selchans_all{10},:); mean_conds_all{11}(selchans_all{11},:); mean_conds_all{12}(selchans_all{12},:)];
plot_mask=[];
for i=ppoi
    plot_mask=cat(1,plot_mask,repmat(reshape(mean_amps_all{i},[size(mean_amps_all{i},1) 1 6]),[1 length(t_ax_short) 1]));
end
plot_mask(find(~isnan(plot_mask)))=1;
plot_data=reshape(plot_data,[size(plot_data,1) 45 6]);
plot_data=plot_data.*plot_mask;
for i=1:size(plot_data,2)
    to_anova=squeeze(plot_data(:,i,:));
    anova_rate=[ones(size(to_anova,1),2) ones(size(to_anova,1),2)*2 ones(size(to_anova,1),2)*3];
    anova_stim=[ones(size(to_anova,1),1) ones(size(to_anova,1),1)*2 ones(size(to_anova,1),1) ones(size(to_anova,1),1)*2 ones(size(to_anova,1),1) ones(size(to_anova,1),1)*2];
    anova_rats=[ones(length(selchans_all{5}),6); ones(length(selchans_all{8}),6)*2; ones(length(selchans_all{9}),6)*3; ones(length(selchans_all{10}),6)*4; ones(length(selchans_all{11}),6)*5; ones(length(selchans_all{12}),6)*6];
    anova_chan=[repmat(1:length(selchans_all{5}),[6 1])'; repmat(1:length(selchans_all{8}),[6 1])'; repmat(1:length(selchans_all{9}),[6 1])'; repmat(1:length(selchans_all{10}),[6 1])'; repmat(1:length(selchans_all{11}),[6 1])'; repmat(1:length(selchans_all{12}),[6 1])'];
    all_ps(i,:)=anovan(to_anova(:),{anova_stim(:) anova_rate(:) anova_rats(:)},'varnames',{'stim' 'rate' 'rat'},'random',[3],'model','full','display','off');
end
all_ps=all_ps(:,[1 2 4]); % stim, rate, stim*rate
for i=1:size(all_ps,2)
    sorted_p=sort(all_ps(:,i));
    p_fdr=sorted_p*0;
    for j=1:length(p_fdr)
        if sorted_p(j)<=.025*j/length(p_fdr)
            p_fdr(j)=1;
        end
    end
    if length(find(p_fdr==1))>0
        p_thresh(i)=sorted_p(max(find(p_fdr==1)));
    else
        p_thresh(i)=0;
    end
end
p_thresh
hold on
shadedErrorBar(t_ax(-1+find(t_ax>=timewin(1)&t_ax<=timewin(2))),nanmean(plot_data(:,:,1),1),nanstd(plot_data(:,:,1),[],1)/sqrt(size(plot_data,1)),{'-','color',[0 0 .5]},1);
shadedErrorBar(t_ax(-1+find(t_ax>=timewin(1)&t_ax<=timewin(2))),nanmean(plot_data(:,:,2),1),nanstd(plot_data(:,:,2),[],1)/sqrt(size(plot_data,1)),{'--','color',[0 0 .5]},1);
shadedErrorBar(t_ax(-1+find(t_ax>=timewin(1)&t_ax<=timewin(2))),nanmean(plot_data(:,:,3),1),nanstd(plot_data(:,:,3),[],1)/sqrt(size(plot_data,1)),{'-','color',[0 0 1]},1);
shadedErrorBar(t_ax(-1+find(t_ax>=timewin(1)&t_ax<=timewin(2))),nanmean(plot_data(:,:,4),1),nanstd(plot_data(:,:,4),[],1)/sqrt(size(plot_data,1)),{'--','color',[0 0 1]},1);
shadedErrorBar(t_ax(-1+find(t_ax>=timewin(1)&t_ax<=timewin(2))),nanmean(plot_data(:,:,5),1),nanstd(plot_data(:,:,5),[],1)/sqrt(size(plot_data,1)),{'-','color',[.5 .5 1]},1);
shadedErrorBar(t_ax(-1+find(t_ax>=timewin(1)&t_ax<=timewin(2))),nanmean(plot_data(:,:,6),1),nanstd(plot_data(:,:,6),[],1)/sqrt(size(plot_data,1)),{'--','color',[.5 .5 1]},1);
xlim([timewin(1)+.01 timewin(2)-.01])
ylabel('amplitude')
xlabel('time (s)')
line([-.2 1],[0 0],'linestyle','--','color',[.5 .5 .5])
line([0 0],[-.2 .8],'linestyle','--','color',[.5 .5 .5])