https://github.com/ViCCo-Group/THINGS-data
Raw File
Tip revision: 2d95c15d3a2cc5984ffd4a9a2c4ad3496847ca9d authored by Oliver Contier on 28 February 2023, 15:15:53 UTC
fixed howtocite
Tip revision: 2d95c15
step2b_data_quality-ERFs.py
#!/usr/bin/env python3

"""
@ Lina Teichmann

    INPUTS: 
    call from command line with following inputs: 
        -bids_dir

    OUTPUTS:
    Plots the ERFs of the repeat trials.

    NOTES:
    If it doesn't exist, the script makes a figures folder in the BIDS derivatives folder
  
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import mne, os
import pandas as pd 
import seaborn as sns


#*****************************#
### PARAMETERS ###
#*****************************#
n_participants              = 4
n_sessions                  = 12
n_images                    = 200
channel_picks               = ['O','T','P']
title_names                 = ['Occipital','Temporal','Parietal']
colors                      = ['mediumseagreen','steelblue','goldenrod','indianred','grey']
plt.rcParams['font.size']   = '16'
plt.rcParams['font.family'] = 'Helvetica'

#*****************************#
### HELPER FUNCTIONS ###
#*****************************#
def load_epochs(preproc_dir,all_epochs = []):
    for p in range(1,n_participants+1):
        epochs = mne.read_epochs(f'{preproc_dir}/preprocessed_P{str(p)}-epo.fif', preload=False)
        all_epochs.append(epochs)
    return all_epochs

# helper function 
def plot_erfs(epochs,n_sessions,name,color,ax,ax2,lab):
    ctf_layout = mne.find_layout(epochs.info)
    picks_epochs = [epochs.ch_names[i] for i in np.where([s[2]==name for s in epochs.ch_names])[0]]
    picks = np.where([i[2]==name for i in ctf_layout.names])[0]

    # get evoked data
    for s in range(n_sessions):    
        evoked = epochs[(epochs.metadata['trial_type']=='test') & (epochs.metadata['session_nr']==s+1)].average()
        evoked.pick_channels(ch_names=picks_epochs)
        ax.plot(epochs.times*1000,np.mean(evoked.data.T,axis=1),color=color,lw=0.5,alpha=0.4)
    evoked = epochs[(epochs.metadata['trial_type']=='test')].average()
    evoked.pick_channels(ch_names=picks_epochs)

    # plot ERFs for selected sensor group
    ax.plot(epochs.times*1000,np.mean(evoked.data.T,axis=1),color=color,lw=1,label=lab)
    ax.set_xlim([epochs.times[0]*1000,epochs.times[len(epochs.times)-1]*1000])
    ax.set_ylim([-0.6,0.6])
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)

    #  plot sensor locations
    ax2.plot(ctf_layout.pos[:,0],ctf_layout.pos[:,1],color='gainsboro',marker='.',linestyle='',markersize=5)
    ax2.plot(ctf_layout.pos[picks,0],ctf_layout.pos[picks,1],color='grey',marker='.',linestyle='',markersize=5)
    ax2.set_aspect('equal')
    plt.axis('off')

#  Make the ERF plot
def make_figure(all_epochs,fig_dir):
    fig = plt.figure(num=1,tight_layout=True,figsize = (11,6))
    gs = GridSpec(3, 5, figure=fig)
    for i,ch in enumerate(channel_picks):
        for p in range(n_participants):
            ax = fig.add_subplot(gs[i, p])
            if i == 0:
                ax.set_title('M' + str(p+1))
            if i == 2:
                ax.set_xlabel('time (ms)')
            else: 
                plt.setp(ax.get_xticklabels(), visible=False)
            if p == 0:
                ax.set_ylabel('fT')
            else: 
                plt.setp(ax.get_yticklabels(), visible=False)

            ax2=fig.add_subplot(gs[i, -1])
        
            plot_erfs(all_epochs[p],12,ch,colors[p],ax,ax2,'Sub' + str(p+1))
        ax2.set_title(title_names[i])
    plt.savefig(f'{fig_dir}/data_quality-ERFs.pdf',dpi=1000)


#*****************************#
### COMMAND LINE INPUTS ###
#*****************************#
if __name__=='__main__':
    import argparse
    parser = argparse.ArgumentParser()

    parser.add_argument(
        "-bids_dir",
        required=True,
        help='path to bids root',
    )

    args = parser.parse_args()
    bids_dir                    = args.bids_dir
    preproc_dir                 = f'{bids_dir}/derivatives/preprocessed/'
    sourcedata_dir              = f'{bids_dir}/sourcedata/'
    fig_dir                      = f'{bids_dir}/derivatives/figures/'
    if not os.path.exists(fig_dir):
        os.makedirs(fig_dir)

    ####### Run ########
    all_epochs = load_epochs(preproc_dir)
    make_figure(all_epochs,fig_dir)
back to top