Browse Source

gin commit from chris-mh16

New files: 1
Modified files: 5
Chris Klink 1 year ago
parent
commit
39029dca63

+ 151 - 25
CODE/figfunctions/figure_SRT.m

@@ -16,7 +16,11 @@ maxRT = 750;
 fastRT = 250;
 
 frt = figure;
-set(frt,'Position',[1200,100,1650,800]);
+set(frt,'Position',[0,0,1900,800]);
+frt2=figure;
+set(frt2,'Position',[0,0,1800,700]);
+
+
 for mi=1:2
     fastRT = ceil(prctile([RTcoll{mi,1};RTcoll{mi,3}],25));   
     fastRT = 10*ceil(fastRT/10);
@@ -73,10 +77,67 @@ for mi=1:2
     xlabel('RT (ms)');ylabel('Probability'); 
     title(['M' num2str(mi) ': fast sacc ' num2str(minRT) '-' num2str(fastRT) ' ms'])
 
+    figure(frt2)
+    subplot(2,4,4+mi);
+    hold on;
+    for c = [3 1]
+        rts = RTcoll{mi,c};
+        rts(rts>fastRT) = [];
+        be=minRT:5:fastRT;
+        histogram(rts,be,'FaceColor',cols(c,:),'Normalization',...
+            'probability'); hold all
+    end
+    set(gca,'xlim',[130 fastRT],'ylim',[0 0.29],'Box','off');
+    xlabel('RT (ms)');ylabel('Probability'); 
+    title(['M' num2str(mi) ': fast sacc ' num2str(minRT) '-' num2str(fastRT) ' ms'])
+
+    subplot(2,4,mi);
+    hold on;
+    nSRT = numel(RTcoll{mi,1})+numel(RTcoll{mi,2})+numel(RTcoll{mi,3});
+    for c = [3 1]
+        rts = RTcoll{mi,c};
+        be=minRT:10:710;
+        histogram(rts,be,'FaceColor',cols(c,:),'Normalization',...
+            'probability'); hold all
+    end
+    %set(gca,'xlim',[130 715],'ylim',[0 0.29],'Box','off');
+    xlabel('RT (ms)');ylabel('Probability'); 
+    title(['M' num2str(mi) ': all sacc (norm by choice)'])
+
+    subplot(2,4,2+mi);
+    hold on;
+    nSRT = numel(RTcoll{mi,1})+numel(RTcoll{mi,2})+numel(RTcoll{mi,3});
+    for c = [3 1]
+        rts = RTcoll{mi,c};
+        be=minRT:10:710;
+        H = histogram(rts,be,'FaceColor',cols(c,:),'Normalization',...
+            'count'); hold all
+        set(H,"BinCounts",H.BinCounts./nSRT);
+    end
+    %set(gca,'xlim',[130 715],'ylim',[0 0.29],'Box','off');
+    xlabel('RT (ms)');ylabel('Probability'); 
+    title(['M' num2str(mi) ': all sacc (norm by all)'])
+
+    subplot(2,4,6+mi);
+    hold on;
+    for c = [3 1]
+        rts = RTcoll{mi,c};
+        rts(rts>fastRT) = [];
+        be=minRT:5:fastRT;
+        H2 = histogram(rts,be,'FaceColor',cols(c,:),'Normalization',...
+            'count'); hold all
+        set(H2,"BinCounts",H2.BinCounts./nSRT);
+    end
+    set(gca,'xlim',[130 fastRT],'ylim',[0 0.13],'Box','off');
+    xlabel('RT (ms)');ylabel('Probability'); 
+    title(['M' num2str(mi) ': fast sacc ' num2str(minRT) '-' num2str(fastRT) ' ms'])
+
+
+
     %STATS
     % choices
     fprintf('\n==============================\n')
-    fprintf(['Stats (25% fastest RTs) for M' num2str(mi) '\n']);
+    fprintf(['Stats (25prc fastest RTs) for M' num2str(mi) '\n']);
     fprintf('==============================\n')
     rt1 = RTcoll{mi,1}; rt1=rt1(rt1<fastRT & rt1>minRT);
     rt2 = RTcoll{mi,2}; rt2=rt2(rt2<fastRT & rt2>minRT);
@@ -114,26 +175,33 @@ for mi=1:2
 end
 
 %% RT_SlidingWindow ====
-aRT = [RTcoll{1,1} ones(size(RTcoll{1,1})) ones(size(RTcoll{1,1}));...
+aRT = [...
+    RTcoll{1,1} ones(size(RTcoll{1,1})) ones(size(RTcoll{1,1}));...
     RTcoll{2,1} 2.*ones(size(RTcoll{2,1})) ones(size(RTcoll{2,1}));...
+    RTcoll{1,2} ones(size(RTcoll{1,2})) 2.*ones(size(RTcoll{1,2}));...
+    RTcoll{2,2} 2.*ones(size(RTcoll{2,2})) 2.*ones(size(RTcoll{2,2}));...
     RTcoll{1,3} ones(size(RTcoll{1,3})) 3.*ones(size(RTcoll{1,3}));...
     RTcoll{2,3} 2.*ones(size(RTcoll{2,3})) 3.*ones(size(RTcoll{2,3}))];
 
-fsw=figure;
-set(fsw,'Position',[1200,100,600,400])
 DoPooled = false;
+fsw=figure;
+set(fsw,'Position',[0,0,900,350]);
+if DoPooled
+    set(fsw,'Position',[0,0,1350,350]);
+end
 
 ws = 20; st = 10; 
 minw = 0; maxw = 500;
 
-xmin=120; %minw+st/2; 
+xmin=125; %minw+st/2; 
 xmax=305; %455;
 
 cw = [minw minw+ws];
-dotsz = 80;
+dotsz = 90;
 
 cRT_T = aRT(aRT(:,3)==1,1);
 cRT_SD = aRT(aRT(:,3)==3,1);
+cRT_ND = aRT(aRT(:,3)==2,1);
 
 medRT = median(aRT(:,1));
 p25 = prctile(aRT(:,1),25);
@@ -142,14 +210,25 @@ mRT = mean(aRT(:,1));
 
 smoothwin = 3;
 
-pSD = [];
+pSD = []; pND = []; pSDR = [];
 while cw(2) < maxw
     TSDratio = ...
         sum(cRT_SD>cw(1) & cRT_SD<=cw(2))./...
         (sum(cRT_T>cw(1) & cRT_T<=cw(2)) + sum(cRT_SD>cw(1) & cRT_SD<=cw(2)));
-    if ~isnan(TSDratio)
+    TNDratio = ...
+        sum(cRT_ND>cw(1) & cRT_ND<=cw(2))./...
+        (sum(cRT_T>cw(1) & cRT_T<=cw(2)) + sum(cRT_ND>cw(1) & cRT_ND<=cw(2)));
+    SDR = sum(cRT_SD>cw(1) & cRT_SD<=cw(2))./size(aRT,1);
+    
+    if ~isnan(TSDratio) 
         pSD=[pSD; mean(cw) TSDratio];
     end
+    if ~isnan(TNDratio)
+        pND=[pND; mean(cw) TNDratio];
+    end
+    if ~isnan(SDR)
+        pSDR=[pSDR; mean(cw) SDR];
+    end
     cw=cw+st;
 end
 if DoPooled
@@ -158,12 +237,7 @@ if DoPooled
     plot([p25 p25],[-0.1 1.1],'k--');
     plot([medRT medRT],[-0.1 1.1],'k-');
     plot([p75 p75],[-0.1 1.1],'k--');
-    
-    x = pSD(pSD(:,1)>xmin & pSD(:,1)<=xmax ,1);
-    y = pSD(pSD(:,1)>xmin & pSD(:,1)<=xmax ,2);
-    xx = x(1):2:x(end);
-    yy = smooth(spline(x,y,xx),smoothwin);
-    scatter(pSD(:,1),pSD(:,2),'ko',...
+    scatter(pSDR(:,1),pSDR(:,2),'ko',...
         'MarkerFaceColor','k','MarkerFaceAlpha',0.9,'SizeData',dotsz);
     set(gca,'xlim',[xmin xmax],'ylim',[-0.02 1.02],'Box','off')
     title('Pooled');
@@ -175,19 +249,31 @@ for mi=1:2
 
     cRT_T = aRT(aRT(:,2)==mi & aRT(:,3)==1,1); 
     cRT_SD = aRT(aRT(:,2)==mi & aRT(:,3)==3,1); 
+    cRT_ND = aRT(aRT(:,2)==mi & aRT(:,3)==2,1); 
+    cRT_A = aRT(aRT(:,2)==mi ,1); 
 
     medRT = median(aRT(aRT(:,2)==mi,1));
     p25 = prctile(aRT(aRT(:,2)==mi,1),25);
     p75 = prctile(aRT(aRT(:,2)==mi,1),75);
     mRT = mean(aRT(aRT(:,2)==mi,1));
 
-    pSD = [];
+    pSDT = [];pSDND = []; pSDR = [];
     while cw(2) < maxw
-        TSDratio = ...
+        SDTratio = ...
             sum(cRT_SD>cw(1) & cRT_SD<=cw(2))./...
             (sum(cRT_T>cw(1) & cRT_T<=cw(2)) + sum(cRT_SD>cw(1) & cRT_SD<=cw(2)));
-        if ~isnan(TSDratio)
-            pSD=[pSD; mean(cw) TSDratio];
+        SDNDratio = ...
+            sum(cRT_SD>cw(1) & cRT_SD<=cw(2))./...
+            (sum(cRT_SD>cw(1) & cRT_SD<=cw(2)) + sum(cRT_ND>cw(1) & cRT_ND<=cw(2))/4);
+        SDR = sum(cRT_SD>cw(1) & cRT_SD<=cw(2))./(sum(cRT_A>cw(1) & cRT_A<=cw(2)));
+        if ~isnan(SDTratio)
+            pSDT=[pSDT; mean(cw) SDTratio];
+        end
+        if ~isnan(SDNDratio)
+            pSDND=[pSDND; mean(cw) SDNDratio];
+        end
+        if ~isnan(SDR)
+            pSDR=[pSDR; mean(cw) SDR];
         end
         cw=cw+st;
     end
@@ -200,13 +286,53 @@ for mi=1:2
     plot([p25 p25],[-0.1 1.1],'k--');
     plot([medRT medRT],[-0.1 1.1],'k-');
     plot([p75 p75],[-0.1 1.1],'k--');
-    x = pSD(pSD(:,1)>xmin & pSD(:,1)<=xmax ,1);
-    y = pSD(pSD(:,1)>xmin & pSD(:,1)<=xmax ,2);
-    xx = x(1):2:x(end);
-    yy = smooth(spline(x,y,xx),smoothwin);
-    scatter(pSD(:,1),pSD(:,2),'ko',...
-        'MarkerFaceColor','k','MarkerFaceAlpha',0.9,'SizeData',dotsz);
+    scatter(pSDR(:,1),pSDR(:,2),'ko',...
+        'MarkerFaceColor','k','MarkerFaceAlpha',0.8,'SizeData',dotsz);
     set(gca,'xlim',[xmin xmax],'ylim',[-0.02 1.02],'Box','off')
     title(['M' num2str(mi)])
     xlabel('RT (ms)'); ylabel('SD choice-ratio')
+    legend({'','p25','median','p75','SD/ALL'});
+end
+
+%% stats on SD ratio
+% proportions SD choices per octile
+nbins=8; 
+for m=1:2
+    mo(m).cnt=[];
+    for i=1:nbins
+        STDR(m,i).m = m;
+        STDR(m,i).prct = i;
+        STDR(m,i).prctval =  prctile(aRT(aRT(:,2)==m,1),i*100/nbins);
+        if i == 1
+            STDR(m,i).SRT_SD = aRT( aRT(:,3)==3 & aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval,1);
+            STDR(m,i).SRT_A = aRT( aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval,1);
+            STDR(m,i).SRT_SDT = aRT( aRT(:,3)~=2 & aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval,1);
+        else
+            STDR(m,i).SRT_SD = aRT(aRT(:,3)==3 & aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval & aRT(:,1)> STDR(m,i-1).prctval,1);
+            STDR(m,i).SRT_A = aRT(aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval & aRT(:,1)> STDR(m,i-1).prctval,1);
+            STDR(m,i).SRT_SDT = aRT(aRT(:,3)~=2 & aRT(:,2)==m & aRT(:,1)< STDR(m,i).prctval & aRT(:,1)> STDR(m,i-1).prctval,1);
+        end
+        STDR(m,i).cnt = [length(STDR(m,i).SRT_SD) length(STDR(m,i).SRT_A) length(STDR(m,i).SRT_SDT)];
+        mo(m).cnt=[mo(m).cnt; STDR(m,i).cnt];
+    end
+
+    % Observed data
+    n1 = sum(mo(m).cnt(1,1)); % fastest 12.5%
+    N1 = sum(mo(m).cnt(1,2));
+    n2 = sum(mo(m).cnt(2:4,1)); % 12.5 - 25%
+    N2 = sum(mo(m).cnt(2:4,2));
+    % Pooled estimate of proportion
+    p0 = (n1+n2) / (N1+N2);
+    % Expected counts under H0 (null hypothesis)
+    n10 = N1 * p0; n20 = N2 * p0;
+    % Chi-square test, by hand
+    observed = [n1 N1-n1 n2 N2-n2];
+    expected = [n10 N1-n10 n20 N2-n20];
+    chi2stat = sum((observed-expected).^2 ./ expected);
+    p = 1 - chi2cdf(chi2stat,1);
+
+    fprintf(['== Monkey ' num2str(m) ' pSD ===\n'])
+    fprintf(['FASTEST 1/8 vs 1/8 to median\n'])
+    fprintf(['CHISQ (1,' num2str(mo(m).cnt(1,2)) ') = ' ...
+        num2str(chi2stat) ', p = ' num2str(p) '\n'])
 end

+ 16 - 0
CODE/figfunctions/figure_binmod.m

@@ -104,6 +104,10 @@ for mi = 1:2
     title([monkeys{mi} ' TARG']); xlabel('time(ms)'); ylabel('dMUA');
     set(gca,'ylim',[-0.04 0.14]);
 
+    firstsig = find(p<0.05/bonf & mn>0,1,'first');
+    fprintf(['M' num2str(mi) '-T: First sign. positive bin at ' ...
+        num2str(bb_c(firstsig)) '\n'])
+
     subplot(2,3,mi+4); hold on;
     df = trcs{3} - trcs{2}; % difference for each channel separately
     mn = mean(df);
@@ -123,6 +127,10 @@ for mi = 1:2
     errorbar(bb_c,mn,csem,'LineStyle','none','Color','k');
     title([monkeys{mi} ' SD']); xlabel('time(ms)'); ylabel('dMUA');
     set(gca,'ylim',[-0.12 0.05]);
+
+    firstsig = find(p<0.05/bonf & mn>0,1,'first');
+    fprintf(['M' num2str(mi) '-SD: First sign. positive bin at ' ...
+        num2str(bb_c(firstsig)) '\n'])
 end
 
 % pooled animals
@@ -158,6 +166,10 @@ errorbar(bb_c,mn,csem,'LineStyle','none','Color','k');
 title('BOTH monkeys TARG'); xlabel('time(ms)'); ylabel('dMUA');
 set(gca,'ylim',[-0.04 0.14]);
 
+firstsig = find(p<0.05/bonf & mn>0,1,'first');
+fprintf(['POOLED-T: First sign. positive bin at ' ...
+    num2str(bb_c(firstsig)) '\n'])
+
 subplot(2,3,4);hold on
 df = trcs{3} - trcs{2}; % difference for each channel separately
 mn = mean(df);
@@ -179,6 +191,10 @@ title('BOTH monkeys SD'); xlabel('time(ms)'); ylabel('dMUA');
 set(gca,'ylim',[-0.12 0.05]);
 text(gca,10,-0.1,'Colored bars: p<0.05, Bonf. corr.')
 
+firstsig = find(p<0.05/bonf & mn>0,1,'first');
+fprintf(['POOLED-SD: First sign. positive bin at ' ...
+    num2str(bb_c(firstsig)) '\n'])
+
 %% save figure
 figure(f1)
 saveas(gcf,fullfile(savedir,'figure_binmod.svg'));

+ 479 - 0
CODE/figfunctions/figure_neurotrace.asv

@@ -0,0 +1,479 @@
+function figure_neurotrace(fld) 
+
+fprintf('\n======================================================\n');
+fprintf('-- Creating Figure neurotrace --\n');
+fprintf('======================================================\n');
+
+%% settings
+green = [23, 105, 13]./255;
+red = [201, 0, 34]./255;
+cols = [green; 0.3 0.3 0.3; red; 1.00,0.84,0.00];
+
+snrthres = 2.5;
+smoothfact = 15; mindays = 3;
+falpha = 1; % set to 1 for converting to illustrator
+
+do_cleantraces = true; % if true we throw out artefact traces
+% citerion: avg emplitude in first 100ms after stimulus onset should be 
+% at least 3 times the standard deviation in the preceding 100 ms
+do_latency = true;
+
+%% load RTs
+if exist(fullfile(fld.basedir, 'results','figure_behavior','RTs.mat'),'file')
+    load(fullfile(fld.basedir, 'results','figure_behavior','RTs.mat'),'RTs');  
+else
+    error('RTs missing, run figure_behavior.m first');
+end
+
+%% plots
+% load traces
+monkeys = {'M1','M2'};
+wm = [];
+
+%% plot
+f1 = figure;
+set(f1,'Position',[500 500 1600 800])
+% individual animals
+for mi = 1:2
+    monkey = monkeys{mi};
+    savedir = fullfile(fld.basedir,'results','figure_neurotrace');
+    load(fullfile(savedir, [monkey '_averages_snr' num2str(snrthres) '_mindays' num2str(mindays) '.mat']),...
+        'traces','tracesLUT','grandtraces','grandtracesLUT','env_t');
+
+    traces_ALL{mi} = traces; %#ok<*AGROW,*NASGU> 
+    tracesLUT_ALL{mi} = tracesLUT;
+    grandtraces_ALL{mi} = grandtraces;
+    grandtracesLUT_ALL{mi} = grandtracesLUT;
+    env_t_ALL{mi} = env_t;
+    
+    figure(f1); subplot(2,3,mi+1); hold on;
+    mns = []; sems = []; trcs = cell(0);
+    for cond = 1:3
+                    
+        if do_cleantraces
+            incl = grandtracesLUT(:,2)==cond;
+            get_traces = grandtraces(incl,:);
+            sel_LUT = grandtracesLUT(incl,:);
+            incl_ALL{mi,cond} = incl;
+
+            prewin = env_t > -0.1 & env_t < 0;
+            postwin = env_t > 0 & env_t < 0.1;
+            sd_pre = std(get_traces(:,prewin),0,2);
+            m_post = mean(get_traces(:,postwin),2);
+            snr_sel = m_post > 3*sd_pre;
+
+            incl = incl(snr_sel');
+            get_traces = get_traces(snr_sel',:);
+            newLUT = sel_LUT(snr_sel',:);
+
+            % now average sites over sessions
+            get_traces_ = get_traces; % backup
+            get_traces = [];
+            for rs = unique(newLUT(:,1))'
+                if sum(newLUT(:,1)==rs)>1
+                    get_traces = [get_traces;...
+                        mean(get_traces_(newLUT(:,1)==rs,:))];
+                else
+                    get_traces = [get_traces; get_traces_(newLUT(:,1)==rs,:)];
+                end
+            end
+            recsites{mi,cond}=unique(newLUT(:,1));
+        else 
+            incl = tracesLUT(:,2)==cond;
+            get_traces = traces(incl,:);
+            sel_LUT = tracesLUT(incl,:);
+            incl_ALL{mi,cond} = incl;
+
+            newLUT = sel_LUT;
+            recsites{mi,cond}=unique(newLUT(:,1));
+        end
+
+        % calculate the mean across channels        
+        s_traces = [];
+        for t = 1:size(get_traces,1)
+            s_traces(t,:) = smooth(get_traces(t,:),smoothfact);
+        end
+        s_traces_ALL{mi,cond} = s_traces;
+
+        mn = mean(s_traces);
+        mns(cond,:) = mn; % for calculating the difference between conditions
+        trcs{cond} = s_traces; % for calculating the modulation onset
+        sem = std(s_traces)./sqrt(sum(incl));
+        sems(cond,:) = sem;
+
+        t2 = [env_t*1e3, fliplr(env_t*1e3)];
+        inBetween = [mn+sem, fliplr(mn-sem)];
+        fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
+            'FaceAlpha',falpha); hold all
+        plot(env_t*1e3,mn,'Color',cols(cond,:));
+        
+        RT=RTs{mi}(cond);
+        plot([RT RT],[0 1],'--','Color',cols(cond,:),'LineWidth',2)
+        
+        xlim([0 0.25*1e3]); ylim([0 1]);
+    end
+    title(monkeys{mi}); xlabel('time(ms)'); ylabel('MUA'); 
+
+
+    % only include sites with all conditions 
+    if size(recsites{mi,1},1) <= size(recsites{mi,2},1) && ...
+            size(recsites{mi,1},1) <= size(recsites{mi,3},1)
+        rsidx = [1 2 3];
+    elseif size(recsites{mi,2},1) < size(recsites{mi,1},1) && ...
+            size(recsites{mi,2},1) < size(recsites{mi,3},1)
+        rsidx = [2 1 3];
+    elseif size(recsites{mi,3},1) < size(recsites{mi,1},1) && ...
+            size(recsites{mi,3},1) < size(recsites{mi,2},1)
+        rsidx = [3 1 2];
+    end
+    rsinc{mi,rsidx(1)} = ismember(recsites{mi,rsidx(1)},recsites{mi,rsidx(1)});
+    rsinc{mi,rsidx(2)} = ismember(recsites{mi,rsidx(2)},recsites{mi,rsidx(1)});
+    rsinc{mi,rsidx(3)} = ismember(recsites{mi,rsidx(3)},recsites{mi,rsidx(1)});
+        
+    inRF = {'target','','distractor'};
+    for cond = [1 3 4]
+        figure(f1)
+        subplot(2,3,mi+4);
+        if cond < 4
+            df = trcs{cond}(rsinc{mi,cond},:) - ...
+                trcs{2}(rsinc{mi,2},:); % difference for each channel separately
+            mn = mean(df);
+            csem = std(trcs{cond}(rsinc{mi,cond},:) - ...
+                trcs{2}(rsinc{mi,2},:))./...
+                sqrt(size(trcs{cond}(rsinc{mi,cond},:),1));
+        else
+            df = trcs{1}(rsinc{mi,1},:) - ...
+                trcs{3}(rsinc{mi,3},:);
+            mn = mean(df);
+            csem = std(trcs{1}(rsinc{mi,1},:) - ...
+                trcs{3}(rsinc{mi,3},:))./...
+                sqrt(size(trcs{1}(rsinc{mi,1},:),1));
+        end
+            
+        inBetween = [mn+csem, fliplr(mn-csem)];
+        fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
+            'FaceAlpha',falpha); hold all
+        plot(env_t*1e3,mn,'Color',cols(cond,:));
+        xlim([0 0.25*1e3]); ylim([-0.15 0.25]);
+        
+        % let's calculate the latency of modulation per channel
+        if do_latency
+            do_latfig=0;
+            if mi==1
+                until = find(env_t<0.25,1,'last');
+            else
+                until = find(env_t<0.21,1,'last');
+            end
+            lat = []; coeff = []; rs = []; converge = [];
+            for i = 1:size(df,1) % loop channels
+                cur_resp = df(i,1:until)';
+                if cond==3
+                    % the latency function cannot fit downward modulation
+                    cur_resp = cur_resp*-1;
+                end
+                [lat(i),~] = latencyfit_jp(cur_resp,env_t(1:until)'*1000,do_latfig);
+            end
+
+            % latency on average
+            cur_resp = mn(1:until)';
+            if cond==3
+                % the latency function cannot fit downward modulation
+                cur_resp = cur_resp*-1;
+            end
+            lt = latencyfit_jp(smooth(cur_resp,20),env_t(1:until)'*1000,do_latfig);
+
+            if cond < 4
+                rng('default') % for reproducibility
+                nsamp = 100;
+                btstrp_LT = bootstrp(nsamp,@latencyfit_jp,smooth(cur_resp,20),env_t(1:until)'*1000,0);
+                fprintf([monkey ' ' inRF{cond} ' - Bootstrapped latency [' num2str(nsamp) ' samples] =====\n'])
+                fprintf(['mean: ' num2str(mean(btstrp_LT,'omitnan')) ', std:' num2str(std(btstrp_LT,'omitnan')) '\n'])
+                BLT{mi,cond}=btstrp_LT;
+            end
+
+            if ~isnan(lt) % if lt = nan, latencyfit_jp gives no figure, so we shouldn't make a title either
+                if cond < 4
+                    title([monkey ', fit on average response, ' inRF{cond} ' in RF']);
+                else
+                    title([monkey ', fit on average response modulation']);
+                end
+            end
+
+            if cond < 4
+                disp([monkey ', latency of modulation onset when ' inRF{cond} ' is in RF: ' ...
+                    num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
+                    'ms as estimated on the average response)']);
+            else
+                disp([monkey ', latency of targ-distr modulation: ' ...
+                    num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
+                    'ms as estimated on the average response)']);
+            end
+            %disp(['individual channel estimates: ' num2str(lat)]);
+
+            ModLatency{mi,cond} = lt;
+            figure(f1);
+            plot([lt lt],[-0.2 0.25],'--','Color',cols(cond,:),'LineWidth',2)
+            if isnan(lt)
+                plot([round(mean(lat,'omitnan')) round(mean(lat,'omitnan'))],[-0.2 0.25],...
+                    ':','Color',cols(cond,:),'LineWidth',2)
+            end
+        end
+    end
+    figure(f1)
+    plot([0 .25*1e3],[0 0],'k');
+    title(monkeys{mi}); xlabel('time(ms)'); ylabel('Modulation(MUA)'); 
+
+    % save the extra cleaned traces
+    save(fullfile(savedir,'Traces_SNR_trialbased.mat'),...
+        's_traces_ALL','recsites','rsinc');
+
+    %% statistics
+    window_means = [];
+    win_idx = find(env_t>0.15 & env_t<0.20);
+
+    for cond = 1:3
+        if do_cleantraces
+            incl = grandtracesLUT(:,2)==cond;
+            get_traces = grandtraces(incl,win_idx); %#ok<*FNDSB> 
+        else
+            incl = tracesLUT(:,2)==cond;
+            get_traces = traces(incl,win_idx);
+        end
+        incl = tracesLUT(:,2)==cond;
+
+        % calculate the average in a window
+        s_traces = [];
+        for t = 1:size(get_traces,1)
+            s_traces(t,:) = smooth(get_traces(t,:),smoothfact);
+        end
+        window_means(:,cond) = mean(s_traces,2);
+    end
+    anova1(window_means);
+    
+    wm = [wm; window_means];
+    wm2{mi} = window_means;
+end
+
+figure(f1); subplot(2,3,1); hold on;
+% pooled animals
+for cond = 1:3
+    traces = [s_traces_ALL{1,cond};s_traces_ALL{2,cond}];
+
+    % calculate the mean across channels
+    % get_traces = traces(incl,:);
+    get_traces = traces;
+    s_traces = [];
+    for t = 1:size(get_traces,1)
+        s_traces(t,:) = smooth(get_traces(t,:),smoothfact);
+    end
+
+    mn = mean(s_traces);
+    mns(cond,:) = mn; % for calculating the difference between conditions
+    trcs{cond} = s_traces; % for calculating the modulation onset
+    sem = std(s_traces)./sqrt(sum(size(s_traces,1)));
+    sems(cond,:) = sem;
+    
+    t2 = [env_t*1e3, fliplr(env_t*1e3)];
+    inBetween = [mn+sem, fliplr(mn-sem)];
+    fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
+        'FaceAlpha',falpha); hold all
+    
+    RT=RTs{3}(cond);
+    plot([RT RT],[0 1],'--','Color',cols(cond,:),'LineWidth',2)
+    plot(env_t*1e3,mn,'Color',cols(cond,:));
+    xlim([0 0.25*1e3]); ylim([0 1]);
+end
+title('BOTH animals pooled'); xlabel('time(ms)'); ylabel('MUA'); 
+
+inRF = {'target','','distractor'};
+
+for cond=1:3
+    rsinc{3,cond} = [rsinc{1,cond};rsinc{2,cond}];
+end
+
+for cond = [1 3 4]
+    figure(f1)
+    subplot(2,3,4);
+
+    if cond < 4
+        df = trcs{cond}(rsinc{3,cond},:) - ...
+            trcs{2}(rsinc{3,2},:); % difference for each channel separately
+        mn = mean(df);
+        csem = std(trcs{cond}(rsinc{3,cond},:) - ...
+            trcs{2}(rsinc{3,2},:))./...
+            sqrt(size(trcs{cond}(rsinc{3,cond},:),1));
+    else
+        df = trcs{1}(rsinc{3,1},:) - ...
+            trcs{3}(rsinc{3,3},:);
+        mn = mean(df);
+        csem = std(trcs{1}(rsinc{3,1},:) - ...
+            trcs{3}(rsinc{3,3},:))./...
+            sqrt(size(trcs{1}(rsinc{3,1},:),1));
+    end
+    
+    inBetween = [mn+csem, fliplr(mn-csem)];
+    fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
+        'FaceAlpha',falpha); hold all
+    plot(env_t*1e3,mn,'Color',cols(cond,:));
+    xlim([0 0.25*1e3]); ylim([-0.15 0.25]);
+    
+    % let's calculate the latency of modulation per channel
+    if do_latency
+        do_latfig=0;
+        if mi==1
+            until = find(env_t<0.25,1,'last');
+        else
+            until = find(env_t<0.21,1,'last');
+        end
+        lat = []; coeff = []; rs = []; converge = [];
+        for i = 1:size(df,1) % loop channels
+            cur_resp = df(i,1:until)';
+            if cond==3
+                % the latency function cannot fit downward modulation
+                cur_resp = cur_resp*-1;
+            end
+            [lat(i),~] = latencyfit_jp(cur_resp,env_t(1:until)'*1000,do_latfig);
+        end
+        
+        %collect latencty
+        clat{cond} = lat;
+
+        % latency on average
+        do_latfig=1;
+        cur_resp = mn(1:until)';
+        if cond==3
+            % the latency function cannot fit downward modulation
+            cur_resp = cur_resp*-1;
+        end
+        lt = latencyfit_jp(smooth(cur_resp,20),env_t(1:until)'*1000,do_latfig);
+
+        if cond < 4
+            rng('default') % for reproducibility
+            nsamp = 100;
+            btstrp_LT = bootstrp(nsamp,@latencyfit_jp,smooth(cur_resp,20),env_t(1:until)'*1000,0);
+            fprintf([inRF{cond} ' - Bootstrapped latency [' num2str(nsamp) ' samples] =====\n'])
+            fprintf(['mean: ' num2str(mean(btstrp_LT,'omitnan')) ', std:' num2str(std(btstrp_LT,'omitnan')) '\n'])
+            BLT{3,cond}=btstrp_LT;
+        end
+
+        if ~isnan(lt) % if lt = nan, latencyfit_jp gives no figure, so we shouldn't make a title either
+            if cond < 4
+                title(['BOTH monkeys, fit on average response, ' inRF{cond} ' in RF']);
+            else
+                title('BOTH monkeys, fit on average response modulation');
+            end
+        end
+        if cond < 4
+            disp(['BOTH monkeys, latency of modulation onset when ' inRF{cond} ' is in RF: ' ...
+                num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
+                'ms as estimated on the average response)']);
+        else
+            disp(['BOTH monkeys, latency of targ-distr modulation: ' ...
+                num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
+                'ms as estimated on the average response)']);
+        end
+        %disp(['individual channel estimates: ' num2str(lat)]);
+
+        figure(f1)
+        plot([lt lt],[-0.2 0.25],'--','Color',cols(cond,:),'LineWidth',2)
+        if isnan(lt)
+            plot([round(mean(lat,'omitnan')) round(mean(lat,'omitnan'))],[-0.2 0.25],...
+                ':','Color',cols(cond,:),'LineWidth',2)
+        end
+
+    end
+end
+figure(f1)
+plot([0 .25*1e3],[0 0],'k');
+title('ALL animals pooled'); xlabel('time(ms)'); ylabel('Modulation (MUA)'); 
+
+%% do a bootstrapped latency analysis to get some stats
+rng('default') % for reproducibility
+nsamp = 1000;
+btstrp_LAT = bootstrp(nsamp,@nanmean,[clat{1}' clat{3}']);
+[h,p,ci,stats]=ttest(btstrp_LAT(:,1)-btstrp_LAT(:,2));
+fprintf('== BTSTRP comparison of latencies ===\n');
+fprintf(['Mean difference in latency (T-SD): ' ...
+    num2str(mean(btstrp_LAT(:,1)-btstrp_LAT(:,2))) ', SEM: ' ...
+    num2str(std(btstrp_LAT(:,1)-btstrp_LAT(:,2))./sqrt(size(btstrp_LAT,1))) ...
+    'ms\n']);
+fprintf([num2str(sum((btstrp_LAT(:,1)-btstrp_LAT(:,2))<0)) '/' num2str(nsamp) ...
+    ' T before SD, p = ' num2str(sum((btstrp_LAT(:,1)-btstrp_LAT(:,2))>0)/nsamp) ...
+    ' that H0[T not before SD] is rejected \n']);
+fprintf(['ttest t(' num2str(stats.df) ') = ' num2str(stats.tstat) ...
+    ', p = ' num2str(p) '\n']);
+
+%% use the bootstrapped latency estimates for stats
+fprintf('== Comparison of bootstrapped T-SD latencies ===\n');
+for mi=1:3
+    [h,p,ci,stats]=ttest(BLT{mi,1}-BLT{mi,3});
+    nn = sum(~isnan());
+    if mi<3
+        fprintf(['-- Monkey ' num2str(mi) ' --\n']);
+    else
+        fprintf('-- BOTH Monkeys --\n');
+    end
+    fprintf(['Mean difference in latency (T-SD): ' ...
+        num2str(mean(BLT{mi,1}-BLT{mi,3},'omitnan')) ', SEM: ' ...
+        num2str(std(BLT{mi,1}-BLT{mi,3},'omitnan')./sqrt(sum(~isnan(BLT{mi,1})))) ...
+        'ms\n']);
+    fprintf([num2str(sum((BLT{mi,1}-BLT{mi,3})<0)) '/' num2str(100) ...
+        ' T before SD, p = ' num2str(sum((BLT{mi,1}-BLT{mi,3})>0)/100) ...
+        ' that H0[T not before SD] is rejected \n']);
+    fprintf(['ttest t(' num2str(stats.df) ') = ' num2str(stats.tstat) ...
+        ', p = ' num2str(p) '\n']);
+end
+
+%% statistics
+% do a simple paired t-test for each monkey
+fprintf('=============================\n');
+fprintf(' Paired T-test \n');
+fprintf('=============================\n');
+clear h p stats
+for mi = 1:2
+    mns = wm2{mi};
+    p = []; vp=[]; ci=1;
+    for cond = [1 3]
+        [vh,vp(end+1)]=vartest2(mns(:,cond),mns(:,2));
+        [h,p(end+1),~,stats{mi,ci}] = ttest(mns(:,cond), mns(:,2));
+        ci=ci+1;
+    end
+    [vh,vp(end+1)]=vartest2(mns(:,1),mns(:,3));
+    [h,p(end+1),~,stats{mi,3}] = ttest(mns(:,1), mns(:,3));
+        
+    disp([monkeys{mi} ', target response increased: p=' num2str(p(1)) ...
+        ' t = ' num2str(stats{mi,1}.tstat) ', df = ' num2str(stats{mi,1}.df) ... 
+        ' (unequal var p = ' num2str(vp(1)) ')']);
+    disp([monkeys{mi} ', distractor response decreased: p=' num2str(p(2)) ...
+        ' t = ' num2str(stats{mi,2}.tstat) ', df = ' num2str(stats{mi,2}.df) ... 
+        ' (unequal var p = ' num2str(vp(2)) ')']);
+    disp([monkeys{mi} ', targ vs distractor: p=' num2str(p(3)) ...
+        ' t = ' num2str(stats{mi,3}.tstat) ', df = ' num2str(stats{mi,3}.df) ... 
+        ' (unequal var p = ' num2str(vp(3)) ')']);
+end
+
+%Both together
+mns = [wm2{1}; wm2{2}];
+p = [];vp=[]; ci=1;
+for cond = [1 3]
+    [vh,vp(end+1)]=vartest2(mns(:,cond),mns(:,2));
+    [h,p(end+1),~,stats{3,ci}] = ttest(mns(:,cond), mns(:,2));
+    ci=ci+1;
+end
+[vh,vp(end+1)]=vartest2(mns(:,1),mns(:,3));
+[h,p(end+1),~,stats{3,3}] = ttest(mns(:,1), mns(:,3));
+
+disp(['BOTH, target response increased: p=' num2str(p(1)) ...
+    ' t = ' num2str(stats{3,1}.tstat) ', df = ' num2str(stats{3,1}.df) ... 
+    ' (unequal var p = ' num2str(vp(1)) ')']);
+disp(['BOTH, distractor response decreased: p=' num2str(p(2)) ...
+    ' t = ' num2str(stats{3,2}.tstat) ', df = ' num2str(stats{3,2}.df) ... 
+    ' (unequal var p = ' num2str(vp(2)) ')']);
+disp(['BOTH, targ vs distractor: p=' num2str(p(3)) ...
+    ' t = ' num2str(stats{3,3}.tstat) ', df = ' num2str(stats{3,3}.df) ... 
+    ' (unequal var p = ' num2str(vp(3)) ')']);
+
+
+%% save figure
+figure(f1)
+saveas(gcf,fullfile(savedir,'figure_neurotrace.svg'));

+ 41 - 0
CODE/figfunctions/figure_neurotrace.m

@@ -181,6 +181,16 @@ for mi = 1:2
                 cur_resp = cur_resp*-1;
             end
             lt = latencyfit_jp(smooth(cur_resp,20),env_t(1:until)'*1000,do_latfig);
+
+            if cond < 4
+                rng('default') % for reproducibility
+                nsamp = 100;
+                btstrp_LT = bootstrp(nsamp,@latencyfit_jp,smooth(cur_resp,20),env_t(1:until)'*1000,0);
+                fprintf([monkey ' ' inRF{cond} ' - Bootstrapped latency [' num2str(nsamp) ' samples] =====\n'])
+                fprintf(['mean: ' num2str(mean(btstrp_LT,'omitnan')) ', std:' num2str(std(btstrp_LT,'omitnan')) '\n'])
+                BLT{mi,cond}=btstrp_LT;
+            end
+
             if ~isnan(lt) % if lt = nan, latencyfit_jp gives no figure, so we shouldn't make a title either
                 if cond < 4
                     title([monkey ', fit on average response, ' inRF{cond} ' in RF']);
@@ -336,6 +346,16 @@ for cond = [1 3 4]
             cur_resp = cur_resp*-1;
         end
         lt = latencyfit_jp(smooth(cur_resp,20),env_t(1:until)'*1000,do_latfig);
+
+        if cond < 4
+            rng('default') % for reproducibility
+            nsamp = 100;
+            btstrp_LT = bootstrp(nsamp,@latencyfit_jp,smooth(cur_resp,20),env_t(1:until)'*1000,0);
+            fprintf([inRF{cond} ' - Bootstrapped latency [' num2str(nsamp) ' samples] =====\n'])
+            fprintf(['mean: ' num2str(mean(btstrp_LT,'omitnan')) ', std:' num2str(std(btstrp_LT,'omitnan')) '\n'])
+            BLT{3,cond}=btstrp_LT;
+        end
+
         if ~isnan(lt) % if lt = nan, latencyfit_jp gives no figure, so we shouldn't make a title either
             if cond < 4
                 title(['BOTH monkeys, fit on average response, ' inRF{cond} ' in RF']);
@@ -383,6 +403,27 @@ fprintf([num2str(sum((btstrp_LAT(:,1)-btstrp_LAT(:,2))<0)) '/' num2str(nsamp) ..
 fprintf(['ttest t(' num2str(stats.df) ') = ' num2str(stats.tstat) ...
     ', p = ' num2str(p) '\n']);
 
+%% use the bootstrapped latency estimates for stats
+fprintf('== Comparison of bootstrapped T-SD latencies ===\n');
+for mi=1:3
+    [h,p,ci,stats]=ttest(BLT{mi,1}-BLT{mi,3});
+    nn = sum(~isnan(BLT{mi,1}-BLT{mi,3}));
+    if mi<3
+        fprintf(['-- Monkey ' num2str(mi) ' --\n']);
+    else
+        fprintf('-- BOTH Monkeys --\n');
+    end
+    fprintf(['Mean difference in latency (T-SD): ' ...
+        num2str(mean(BLT{mi,1}-BLT{mi,3},'omitnan')) ', SEM: ' ...
+        num2str(std(BLT{mi,1}-BLT{mi,3},'omitnan')./sqrt(nn)) ...
+        'ms\n']);
+    fprintf([num2str(sum((BLT{mi,1}-BLT{mi,3})<0)) '/' num2str(nn) ...
+        ' T before SD, p = ' num2str(sum((BLT{mi,1}-BLT{mi,3})>0)/nn) ...
+        ' that H0[T not before SD] is rejected \n']);
+    fprintf(['ttest t(' num2str(stats.df) ') = ' num2str(stats.tstat) ...
+        ', p = ' num2str(p) '\n']);
+end
+
 %% statistics
 % do a simple paired t-test for each monkey
 fprintf('=============================\n');

+ 6 - 5
CODE/functions/latencyfit_jp.m

@@ -37,7 +37,7 @@ fo_ = fitoptions('method','NonlinearLeastSquares','MaxFunEvals',1000,'Lower',[0,
 set(fo_,'Startpoint',st_);
 
 try
-[cfun,gof,output] = fit(t,y,ft_,fo_);
+    [cfun,gof,output] = fit(t,y,ft_,fo_);
 catch
     %Fit went wrong
     lat = NaN;
@@ -71,7 +71,8 @@ ixs = [1:mfix]; % only search before max peak
 lat = t(l33ix);
 
 if fig
-figure;plot(t,y,'b',t,yf,'r');legend({'data','fit'});
-hold on;plot([t(l33ix),t(l33ix)],get(gca,'YLim')); % point where 33% of max response reached
-end
-
+    figure; hold on;
+    plot(t,y,'b',t,yf,'r');
+    legend({'data','fit'});
+    plot([t(l33ix),t(l33ix)],get(gca,'YLim')); % point where 33% of max response reached
+end

+ 1 - 1
CODE/main.m

@@ -27,7 +27,7 @@ fld.procdatadir = fullfile(fld.basedir,'Processed_Data');
 
 fprintf('=============================================================\n')   
 fprintf(['System used: ' ThisSystem '\n'])
-fprintf(['Basedir used: ' basedir '\n'])      
+fprintf(['Basedir used: ' fld.basedir '\n'])      
 fprintf('=============================================================\n')      
 
 % info --