123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257 |
- 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;
|