NHP-pRF/Analysis_code/cRF/cRF_process.m

%% base folders >> adjust for your system =================================
base_fld = fullfile(pwd, 'cRF');
% NB! Monkey L = M3; Monkey A = M4

%% Do some processing and organization ====================================
% read csv files
warning off
cRF=readtable(fullfile(base_fld,'allRF.csv'));
warning on

minSNR = 3; minR2 = 0.25;

SNR = [...
    cRF.SNR_fromRF_rightward ...
    cRF.SNR_fromRF_leftward ...
    cRF.SNR_fromRF_upward ...
    cRF.SNR_fromRF_downward ];
SNR(SNR<minSNR)=NaN;
cRF.incSNR = mean(SNR,2)>0;

R2 = [...
    cRF.R2_rightward ...
    cRF.R2_leftward ...
    cRF.R2_upward ...
    cRF.R2_downward ];
R2(R2<minR2)=NaN;
cRF.incR2 = mean(R2,2)>0;

cRF.incSNRandR2 = mean(SNR,2)>minSNR & mean(R2,2)>minR2;

% rewrite borders
cRF.Lp = cRF.RF_left_edge_pixels_;
cRF.Rp = cRF.RF_right_edge_pixels_;
cRF.Tp = cRF.RF_top_edge_pixels_;
cRF.Bp = cRF.RF_bottom_edge_pixels_;

% rewrite centers
cRF.HCd = cRF.RFCenterX_degrees_;
cRF.VCd = cRF.RFCenterY_degrees_;

% recalc
cRF.HSp = cRF.Rp-cRF.Lp;
cRF.HCp = (cRF.Rp+cRF.Lp)./2;
cRF.VSp = cRF.Tp-cRF.Bp;
cRF.VCp = (cRF.Tp+cRF.Bp)./2;

% calculate pixel-deg conversion
ppd = nanmean(cRF.HCp./cRF.HCd);

% size in deg
cRF.HSd = cRF.HSp./ppd; % in deg based on 1.65*SD
cRF.VSd = cRF.VSp./ppd; % in deg based on 1.65*SD

% size in SD >> NOTE THAT THIS IS RADIUS NOW, DO NOT /2 LATER!!
cRF.HSsd = cRF.HSd./3.3; % in deg based on SD 
cRF.VSsd = cRF.VSd./3.3; % in deg based on SD 
cRF.SZsd = sqrt((cRF.HSsd.^2)+(cRF.VSsd.^2));

% add area
cRF.Area = ones(size(cRF.Array_ID));
cRF.Area(strcmp(cRF.MID, 'A') & (cRF.Array_ID == 2 | cRF.Array_ID == 5))=4;
cRF.Area(strcmp(cRF.MID, 'L') & (cRF.Array_ID == 2 | cRF.Array_ID == 3))=4;

%% Horizontal/vertical cRF size ===========================================
% horizontal vs vertical size
M3idx = strcmp(cRF.MID,'L') & cRF.incSNRandR2 ;
M4idx = strcmp(cRF.MID,'A') & cRF.incSNRandR2 ;

figure; 
subplot(1,3,1);hold on
plot([-1 6],[-1 6],'r')
scatter(cRF.HSsd(M3idx), cRF.VSsd(M3idx),'b')
scatter(cRF.HSsd(M4idx), cRF.VSsd(M4idx),'k')
plot([0 0],[-1 6],'k--');plot([-1 6],[0 0],'k--');
xlabel('Horizontal RF-size (SD) edge-based')
ylabel('Vertical RF-size (SD) edge-based')
set(gca, 'xlim',[-1 6],'ylim',[-1 6]);

subplot(1,3,2);hold on
plot([-1 6],[-1 6],'r')
scatter(cRF.HSsd(M3idx), cRF.SZsd(M3idx),'b')
scatter(cRF.HSsd(M4idx), cRF.SZsd(M4idx),'k')
plot([0 0],[-1 6],'k--');plot([-1 6],[0 0],'k--');
xlabel('Horizontal RF-size (SD) edge-based')
ylabel('RF-size (SD) edge-based')
set(gca, 'xlim',[-1 6],'ylim',[-.5 6]);


subplot(1,3,3);hold on
plot([-1 6],[-1 6],'r')
scatter(cRF.VSsd(M3idx), cRF.SZsd(M3idx),'b')
scatter(cRF.VSsd(M4idx), cRF.SZsd(M4idx),'k')
plot([0 0],[-1 6],'k--');plot([-1 6],[0 0],'k--');
xlabel('Vertical RF-size (SD) edge-based')
ylabel('RF-size (SD) edge-based')
set(gca, 'xlim',[-1 6],'ylim',[-.5 6]);

legend({'x=y','M3','M4'}, 'Location','SouthEast')
fprintf(['nChan L: ' num2str(sum(M3idx)) ...
    ', nChan A: ' num2str(sum(M4idx)) '\n'])

%% cRF symmetry ===========================================================
figure;
subplot(1,3,1);hold on
plot([-1 6],[-1 6],'k--')
scatter(cRF.HSsd(M3idx), cRF.VSsd(M3idx),'k')
scatter(cRF.HSsd(M4idx), cRF.VSsd(M4idx),'r')
xlabel('Horizontal RF-size (SD)')
ylabel('Vertical RF-size (SD)')
set(gca, 'xlim',[0 5],'ylim',[0 5]);
legend({'HSz=VSz','M3','M4'})
subplot(1,3,2);hold on
histogram(cRF.HSsd(M3idx)./cRF.VSsd(M3idx),0.45:0.05:1.55)
nm=nanmedian(cRF.HSsd(M3idx)./cRF.VSsd(M3idx));
plot([nm nm],[0 120],'k--')
title(['M3, median:' num2str(nm)]); xlabel('HSz/VSz'); ylabel('#electrodes')
subplot(1,3,3);hold on
histogram(cRF.HSsd(M4idx)./cRF.VSsd(M4idx),0.45:0.05:1.55)
nm=nanmedian(cRF.HSsd(M4idx)./cRF.VSsd(M4idx));
plot([nm nm],[0 100],'k--')
title(['M4, median:' num2str(nm)]); xlabel('HSz/VSz'); ylabel('#electrodes')

% stats
[p,h,stats] = signrank(cRF.HSsd(M3idx),cRF.VSsd(M3idx));
fprintf(['M3: Difference HSz vs VSz: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);
[p,h,stats] = signrank(cRF.HSsd(M4idx),cRF.VSsd(M4idx));
fprintf(['M4: Difference HSz vs VSz: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

%% aspect ratio ===========================================================
sigmat3 = sort([cRF.HSsd(M3idx), cRF.VSsd(M3idx)],2,'descend');
sigmat4 = sort([cRF.HSsd(M4idx), cRF.VSsd(M4idx)],2,'descend');

figure;
subplot(1,3,1);hold on
plot([-1 6],[-1 6],'k--')
scatter(sigmat3(:,1), sigmat3(:,2),'k')
scatter(sigmat4(:,1), sigmat4(:,2),'r')
xlabel('Largest SD'); ylabel('Smallest SD')
set(gca, 'xlim',[0 5],'ylim',[0 5]);
legend({'unity','M3','M4'})

subplot(1,3,2);hold on
histogram(sigmat3(:,1)./sigmat3(:,2),0.95:0.05:5.05)
nm=nanmedian(sigmat3(:,1)./sigmat3(:,2));
IQR=[nm-iqr(sigmat3(:,1)./sigmat3(:,2))/2 nm+iqr(sigmat3(:,1)./sigmat3(:,2))/2]; 
plot([nm nm],[0 180],'k--')
title(['M3, median:' num2str(nm)]); xlabel('aspect ratio'); ylabel('#electrodes')
fprintf(['M3 median: ' num2str(nm) ', IQR:' num2str(IQR(1)) ' ' num2str(IQR(2)) '\n'])
ar=sigmat3(:,1)./sigmat3(:,2);
fprintf(['M3: ' num2str(sum(ar>2)) '/' num2str(length(ar)) ...
    ' cRFs with aspect ratio >2 (' num2str(100*sum(ar>2)/length(ar)) '%%)\n' ]);

subplot(1,3,3);hold on
histogram(sigmat4(:,1)./sigmat4(:,2),0.95:0.05:5.05)
nm=nanmedian(sigmat4(:,1)./sigmat4(:,2));
IQR=[nm-iqr(sigmat4(:,1)./sigmat4(:,2))/2 nm+iqr(sigmat4(:,1)./sigmat4(:,2))/2]; 
plot([nm nm],[0 120],'k--')
title(['M4, median:' num2str(nm)]); xlabel('aspect ratio'); ylabel('#electrodes')
fprintf(['M4 median: ' num2str(nm) ', IQR:' num2str(IQR(1)) ' ' num2str(IQR(2)) '\n'])
ar=sigmat4(:,1)./sigmat4(:,2);
fprintf(['M4: ' num2str(sum(ar>2)) '/' num2str(length(ar)) ...
    ' cRFs with aspect ratio >2 (' num2str(100*sum(ar>2)/length(ar)) '%%)\n' ]);

% stats
[p,h,stats] = signrank(sigmat3(:,1), sigmat3(:,2));
fprintf(['M3: Difference HSz vs VSz: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);
[p,h,stats] = signrank(sigmat4(:,1), sigmat4(:,2));
fprintf(['M4: Difference HSz vs VSz: Wilcoxon z = ' ...
    num2str(stats.zval) ', p = ' num2str(p) '\n']);

%% restructure for compatibility
warning off
cRF_map1 = readtable(fullfile(base_fld,'combined_A_RF_with_SD.csv'));
cRF_map2 = readtable(fullfile(base_fld,'combined_L_RF_with_SD.csv'));
cRF_map = [cRF_map1;cRF_map2];
warning on

clear RFs
M = {'M03','M04'}; % L,A

% add NSP number
Arr_num = 1:16;
NSP_num = [1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8];

clear idx
idx{1} = strcmp(cRF.MID,'L');
idx{2} = strcmp(cRF.MID,'A');
    
for m=1:2
    save_fld = fullfile(base_fld,M{m}); 
    mkdir(save_fld)

    Tbl=cRF(idx{m},:);
    
    idxm = strcmp(cRF_map.Monkey,M{m});
    Tbl_map = cRF_map(idxm,:);
    NSP_map=[];
    
    for i = 1:8
        a = Arr_num(NSP_num == i);
                
        arr_sel = (Tbl.Array_ID == a(1) | Tbl.Array_ID == a(2));       

        NSP_data = [Tbl.Electrode_ID(arr_sel) ...
            Tbl.Array_ID(arr_sel) ...
            nan(size(Tbl.Array_ID(arr_sel))) ];
        
        for el = 1:length(NSP_data)
            NSP_data(el,3) = ...
                Tbl_map.within_NSP_electrode_ID(...
                Tbl_map.Electrode_ID == NSP_data(el,1) & ...
                Tbl_map.Array_ID == NSP_data(el,2));
        end
        NSP_order = NSP_data(:,3);
        
        meanChannelSNR = mean([...
            Tbl.SNR_fromRF_rightward(arr_sel) ...
            Tbl.SNR_fromRF_leftward(arr_sel) ...
            Tbl.SNR_fromRF_upward(arr_sel) ...
            Tbl.SNR_fromRF_downward(arr_sel) ...
            ],2);
        meanChannelSNR = meanChannelSNR(NSP_order);
        
        meanChannelR2 = mean([...
            Tbl.R2_rightward(arr_sel) ...
            Tbl.R2_leftward(arr_sel) ...
            Tbl.R2_upward(arr_sel) ...
            Tbl.R2_downward(arr_sel) ...
            ],2);
        meanChannelR2 = meanChannelR2(NSP_order);
        
        XX = Tbl.RFCenterX_degrees_(arr_sel);
        YY = Tbl.RFCenterY_degrees_(arr_sel);
        SZsd = Tbl.SZsd(arr_sel);
        HV = [Tbl.HSsd(arr_sel) Tbl.VSsd(arr_sel)];
        TH = Tbl.RFTheta_degrees_(arr_sel);
        SZlisted =  Tbl.RFSize_degrees_(arr_sel);
        Areas = Tbl.Area(arr_sel);
        
        for c = 1:128
            RFs{1,NSP_order(c)}.centrex = XX(c)*ppd;      
            RFs{1,NSP_order(c)}.centrey = YY(c)*ppd;
            RFs{1,NSP_order(c)}.xdeg = XX(c);      
            RFs{1,NSP_order(c)}.ydeg = YY(c);
            RFs{1,NSP_order(c)}.sz = SZsd(c)*ppd;
            RFs{1,NSP_order(c)}.szdeg = SZsd(c); 
            RFs{1,NSP_order(c)}.szdeg_hv = HV(c,:);
            RFs{1,NSP_order(c)}.szlisteddeg = SZlisted(c); 
            RFs{1,NSP_order(c)}.theta = TH(c);
            RFs{1,NSP_order(c)}.ecc = sqrt(XX(c).^2 + YY(c).^2);
            RFs{1,NSP_order(c)}.area = Areas(c);
        end
        save([save_fld '/RFs_instance' num2str(i)], ...
            'meanChannelSNR','meanChannelR2','RFs')
    end
end