Human_Amygdala_MUA_sEEG_FearVideo/code_MATLAB/Figure_1.m

function fig = Figure_1(file_path, varargin)
% Plot Figure 1
%
% include plotting toolbox gramm: https://github.com/piermorel/gramm
%
% neuron: neurons' number
% SPIKES.spikes: neurons x [time; condition; trial timing]
% SPIKES.spikes_wave: neurons x [mean; std]
%

if(nargin>1)
    nExampleNeuron = varargin{1};
else
    nExampleNeuron = 1;
end

neurons = 1;
for nSubject = 1:9
    file_name = sprintf('Data_Subject_%.2d_Session_01.h5',nSubject);
    f = nix.File([file_path,file_name],nix.FileMode.ReadOnly);
    
    block = f.blocks{1};
    group_SpikeTimes = block.openGroup('Spike times');
    trials_labels = cellfun(@(x) x.name, group_SpikeTimes.dataArrays, ...
        'UniformOutput',0);
    info = cellfun(@(x) strsplit(x,'_'), trials_labels, 'UniformOutput',0);
    units = cell2mat(cellfun(@(x) str2double(x{4}), info, 'UniformOutput',0));
    trials = cell2mat(cellfun(@(x) str2double(x{8}), info, 'UniformOutput',0));
    
    for u=1:max(units)
        trial_ind = [];
        SpikeTimes = [];
        for t=1:max(trials)
            condition = (units==u)&(trials==t);
            dataArray_SpikeTimes = group_SpikeTimes.dataArrays{condition};
            data = dataArray_SpikeTimes.readAllData';
            SpikeTimes = [SpikeTimes; data];
            trial_ind = [trial_ind; ones(length(data), 1)*t];
        end
        SPIKES(neurons).spikes = [SpikeTimes trial_ind];
        indMultiTagUnitSpikeTimes = contains(cellfun(@(x) x.name,block.multiTags, ...
            'UniformOutput',0),['Unit_',num2str(u),'_']);
        multiTag_SpikeTimes = block.multiTags{indMultiTagUnitSpikeTimes};
        dataArray_Waveform = multiTag_SpikeTimes.features{1}.openData;
        SPIKES(neurons).waveform = dataArray_Waveform.readAllData';
        neurons = neurons + 1;
    end
    
end


%%
% FR all
FR_all = [];
for i = 1:length(SPIKES)
    n_trials = max(SPIKES(i).spikes(:, 2));
    FR_all = [ FR_all  length(SPIKES(i).spikes(:,1)) / n_trials*0.026];
end
% ISI
isis = [];
for i = 1:length(SPIKES)
    isis_ = [];
    for t = 1:max(SPIKES(i).spikes(:,2))
        isis_ = [isis_; diff(sort(SPIKES(i).spikes(SPIKES(i).spikes(:,2)==t, 1))*1000)];
    end
    isis = [isis; (length(isis_(isis_<3))/length(isis_))*100];
end

% SNR
snr = [];
waveform = [];
for i = 1:length(SPIKES)
    snr_0 = max(abs(SPIKES(i).waveform(:,1)))/mean(SPIKES(i).waveform(:,2));
    snr = [snr snr_0];
    waveform = [waveform; SPIKES(i).waveform(:,1)'];
end

FR = SPIKES(neuron).spikes;

%% Firing rate one
c = [];
fs = -2; % min(FR(:,3));
ls = 24; % max(FR(:,3));
binSize = 0.1;
for i = 1:max(FR(:,2))
    for t = [0 1] % 1 - red, 0 - blue
        if mod(i,2) == t % odd - neutral, even - aversiv
            spikeTimes = FR(FR(:,2) == i, 1);
            trials{i} = spikeTimes;
            c = [c t];
            temp = histc(spikeTimes, fs:binSize:ls);
            binned{i} = smooth(temp(1:end-1)/binSize, 10);
            % we set here bin centers
            bins{i} = (fs+binSize/2):binSize:(ls-binSize/2);
        end
    end
end

%%
clear g
fig = figure;
% FR
g(1)=gramm( 'x', bins, 'y', binned, 'color', c );
n = length(binned);
custom_statfun = @(y)([mean(y);
    mean(y) - std(y)/sqrt(17);
    mean(y) + std(y)/sqrt(17)]);
g(1).stat_summary('setylim', true, 'type', custom_statfun);
g(1).set_names( 'x', '', 'y', 'Rate [Hz]' );
g(1).set_layout_options('position', [0.05 0.7 0.9 0.3])
g(1).no_legend()
g(1).axe_property('xlim',[fs 20],'XColor','none');
% Inset
g(6)=gramm('x', ([1:64]-20)/32,'y', waveform(neuron,:));
g(6).geom_line();
g(6).set_names( 'x', 'Time (ms)', 'y', 'uV' );
g(6).set_layout_options('position',[0.1 0.9 0.1 0.1])
g(6).set_text_options('base_size',4)
g(6).set_color_options('chroma',0,'lightness',30);
g(6).set_line_options('base_size', 1);
g(6).no_legend()
% g(6).axe_property('xlim',[fs 20],'XColor','none');
% raster
g(2)=gramm('x', trials,'color', c);
g(2).geom_raster( 'geom', 'point');
g(2).set_point_options('base_size', 3);
g(2).set_names( 'x', 'Time [s]', 'y', 'Trial no. (reordered)' );
g(2).set_layout_options('position', [0.05 0.43 0.9 0.3])
g(2).no_legend()
g(2).axe_property('xlim',[fs 20], 'ylim', [-1 9]);
% ISI
g(3)=gramm('x', isis);
g(3).stat_bin();
g(3).set_layout_options('position',[0.05 0.1 0.25 0.3])
g(3).set_names( 'x', 'Percentage of ISI < 3 ms', 'y', 'No. of units' );
g(3).set_color_options('hue_range',[155 200]);
% FR all
g(4)=gramm('x', FR_all);
g(4).stat_bin();
g(4).set_layout_options('position',[0.38 0.1 0.25 0.3])
g(4).set_names( 'x', 'Firing rate [Hz]', 'y', 'No. of units' );
g(4).axe_property('xlim',[0 20]);
g(4).set_color_options('hue_range',[155 200]);
% SNR
g(5)=gramm('x', snr);
g(5).stat_bin();
g(5).set_layout_options('position',[0.71 0.1 0.25 0.3])
g(5).set_names( 'x', 'Waveform peak SNR', 'y', 'No. of neurons' );
g(5).axe_property('xlim',[0 10]);
g(5).set_color_options('hue_range',[155 200]);
% Labels
x = [0 0 0 0.35 0.68];
y = [0.97 0.7 0.37 0.37 0.37];
l = ['a' 'b' 'c' 'd' 'e'];
g(7)=gramm('x', x, 'y', y,'label',l);
g(7).geom_label('Color','k', 'fontweight', 'bold')
g(7).axe_property('XColor','none','YColor','none','xlim',[0 1],...
    'ylim',[0 1],'Visible','off');
g(7).set_layout_options('position',[0 0 1 1])
g.draw();
print('-dtiff','-r300','fig_1.tiff')

end