auksztulewicz_et_al_2023_bmcbiol/scripts/Figure4A.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_LFP.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_lfp{p,j});

        if data_present>0




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

            preceding_all_lfp_bc{p,j}=preceding_all_lfp{p,j}-repmat(nanmean(preceding_all_lfp{p,j}(:,min(find(t_ax>=baseline(1))):max(find(t_ax<=baseline(2))),:),2),[1 size(preceding_all_lfp{p,j},2) 1]);
            omissions_all_lfp_bc{p,j}=omissions_all_lfp{p,j}-repmat(nanmean(omissions_all_lfp{p,j}(:,min(find(t_ax>=baseline(1))):max(find(t_ax<=baseline(2))),:),2),[1 size(omissions_all_lfp{p,j},2) 1]);
            preceding_all_lfp_bc{p,j}=preceding_all_lfp_bc{p,j}./repmat(nanstd(preceding_all_lfp_bc{p,j}(:,min(find(t_ax>=baseline(1))):max(find(t_ax<=baseline(2))),:),[],2),[1 size(preceding_all_lfp_bc{p,j},2) 1]);
            omissions_all_lfp_bc{p,j}=omissions_all_lfp_bc{p,j}./repmat(nanstd(omissions_all_lfp_bc{p,j}(:,min(find(t_ax>=baseline(1))):max(find(t_ax<=baseline(2))),:),[],2),[1 size(omissions_all_lfp_bc{p,j},2) 1]);
        end
    end

    mean_evoked=[];
    for j=1:3
        if length(preceding_all_lfp{p,j})>0
            mean_evoked=[mean_evoked nanmean(preceding_all_lfp_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_lfp{p,j})>0
            preceding_all_lfp_bc{p,j}=preceding_all_lfp_bc{p,j}./repmat(totalrange,[1 size(preceding_all_lfp_bc{p,j},2) size(preceding_all_lfp_bc{p,j},3)]);
            omissions_all_lfp_bc{p,j}=omissions_all_lfp_bc{p,j}./repmat(totalrange,[1 size(omissions_all_lfp_bc{p,j},2) size(omissions_all_lfp_bc{p,j},3)]);
        end
    end

    allstd=[];
    for j=1:size(preceding_all_lfp_bc,2)
        try
            allstd=[allstd; squeeze(mean(std(omissions_all_lfp_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_lfp_bc,2)
        try
            badtrls=find(squeeze(mean(std(omissions_all_lfp_bc{p,j}(setdiff(1:64,badchan),find(t_ax>=timewin(1)&t_ax<=timewin(2)),:),[],2),1))>stdcutoff);
            omissions_all_lfp_bc{p,j}(:,:,badtrls)=NaN;
        catch
        end
    end

    allstd=[];
    for j=1:size(preceding_all_lfp_bc,2)
        try
            allstd=[allstd; squeeze(mean(std(preceding_all_lfp_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_lfp_bc,2)
        try
            badtrls=find(squeeze(mean(std(preceding_all_lfp_bc{p,j}(setdiff(1:64,badchan),find(t_ax>=timewin(1)&t_ax<=timewin(2)),:),[],2),1))>stdcutoff);
            preceding_all_lfp_bc{p,j}(:,:,badtrls)=NaN;
        catch
        end
    end

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


            mean_response=[mean_response nanmean(preceding_all_lfp_bc{p,j}(:,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),:),3) nanmean(omissions_all_lfp_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_lfp_bc{p,j}(:,min(find(t_ax>=timewin(1))):max(find(t_ax<=timewin(2))),:),3),2) size(nanmean(omissions_all_lfp_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_lfp_bc,2)
            try
                alltrls=[alltrls nanmean(squeeze(preceding_all_lfp_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_lfp_bc,2)
            try
                alltrls=[alltrls nanmean(squeeze(omissions_all_lfp_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(find(p_omi<p_fdr_thresh&p_prec<p_fdr_thresh),badchan);
    %     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 9: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

figure;
subplot(2,2,2)
hold on
for i=1:3
    scatter(anova_amps(find(anova_stim==1&anova_rate==i)),anova_amps(find(anova_stim==2&anova_rate==i)),25,colscols(i,:),'filled');
    [B,BINT] = regress(anova_amps(find(anova_stim==2&anova_rate==i)),[anova_amps(find(anova_stim==1&anova_rate==i))*0+1 anova_amps(find(anova_stim==1&anova_rate==i))]);
    mdl=fitlm(anova_amps(find(anova_stim==1&anova_rate==i)),anova_amps(find(anova_stim==2&anova_rate==i)));
    t_cookd = 3*mean(mdl.Diagnostics.CooksDistance,'omitnan');
    outlrs=find(mdl.Diagnostics.CooksDistance > t_cookd);
    outlrs=[];
    Xcorr=anova_amps(find(anova_stim==1&anova_rate==i));
    Ycorr=anova_amps(find(anova_stim==2&anova_rate==i));
    [r,p]=corr(Xcorr(setdiff(1:length(Xcorr),outlrs)),Ycorr(setdiff(1:length(Ycorr),outlrs)))
    if p<.05
        plot([-.5 1.5],[-.5 1.5]*B(2)+B(1),'color',colscols(i,:),'linestyle','-')
    else
        plot([-.5 1.5],[-.5 1.5]*B(2)+B(1),'color',colscols(i,:),'linestyle','--')
    end
end
xlim([-.5 1.5])
ylim([-.5 1.5])
line([-.5 1.5],[0 0],'linestyle','--','color',[.5 .5 .5])
line([0 0],[-.5 1.5],'linestyle','--','color',[.5 .5 .5])
sfh1 = subplot(2,2,4)
hold on
for i=1:3
    a=hist(anova_amps(find(anova_stim==1&anova_rate==i)),-.5:.25:1.5);
    line([-.5:.25:1.5],-a/sum(a),'color',colscols(i,:))
end
xlim([-.5 1.5])
set(gca,'ytick',[])
ylim([-.5 0])
sfh1.Position(2)=sfh1.Position(2)+sfh1.Position(4)/2;
sfh1.Position(4)=sfh1.Position(4)/2;
sfh2 = subplot(2,2,1)
hold on
for i=1:3
    a=hist(anova_amps(find(anova_stim==2&anova_rate==i)),-.5:.25:1.5);
    line(-a/sum(a),[-.5:.25:1.5],'color',colscols(i,:))
end
ylim([-.5 1.5])
set(gca,'xtick',[])
xlim([-.5 0])
sfh2.Position(1)=sfh2.Position(1)+sfh1.Position(3)/2;
sfh2.Position(3)=sfh1.Position(3)/2;