sEEG bipolar w/ MNE-Python

seeg_bip.py

import sys
import glob
import numpy as np
from mne import fiff, event, epochs
from iidc import detect, stat, util, cluster

#### Data

# Load the data, then eliminate all but the sEEG channels

# In[4]:

pdfname = glob.glob(u'/data/nasmeg/DATA/Protocoles_epilepsie/simultane/MEG_SEEG/seeg_006_gou_mi/s_Gou_Mi/*/3/c,rfDC')[0]
data = fiff.bti.read_raw_bti(pdfname)
util.bti_import_fix_channel_names(data, pdfname)
picks = fiff.pick_types(data.info, eeg=True, meg=False, exclude=['VEOG', 'E90', 'E91', 'E92', 'E93', 'E94', 'ECG'])
data.save('goumi-run3-seeg.fif', picks=picks, overwrite=True, verbose=False)
del data

# Out[4]:

#     Reading 4D PDF file /data/nasmeg/DATA/Protocoles_epilepsie/simultane/MEG_SEEG/seeg_006_gou_mi/s_Gou_Mi/02%21%13@1030/3/c,rfDC...
#     Creating Neuromag info structure ...
#     ... Setting channel info structure.
#     ... putting coil transforms in Neuromag coordinates
#     ... Reading digitization points from /data/nasmeg/DATA/Protocoles_epilepsie/simultane/MEG_SEEG/seeg_006_gou_mi/s_Gou_Mi/02%21%13@1030/3/hs_file
#     ... putting digitization points in Neuromag coordinates
#     ... Computing new device to head transform.
#     Done.
#     Warning. Currently direct inclusion of 4D weight tables is not supported. For critical use cases 
#     please take into account the MNE command 'mne_create_comp_data' to include weights as printed out 
#     by the 4D 'print_table' routine.
#     Reading raw data from /data/nasmeg/DATA/Protocoles_epilepsie/simultane/MEG_SEEG/seeg_006_gou_mi/s_Gou_Mi/02%21%13@1030/3/c,rfDC...
#         Range : 0 ... 1453885 =      0.000 ...   714.612 secs
#     Ready.
# 

# In[4]:

seeg = fiff.Raw('goumi-run3-seeg.fif', preload=True)


# Out[4]:

#     Opening raw data file /home/duke/tools/iidc/work/goumi-run3-seeg.fif...
#     Current compensation grade : 0
#         Range : 0 ... 1453885 =      0.000 ...   714.612 secs
#     Ready.
#     Adding average EEG reference projection.
#     Reading 0 ... 1453885  =      0.000 ...   714.612 secs...
#     [done]
# 

# Now, knock off a few that survived the pick, 

# In[12]:

picks = array([i for i, c in enumerate(seeg.info['chs']) if c['kind']==2])
seeg.save('goumi-run3-seeg.fif', picks=picks, overwrite=True, verbose=False)


# In[65]:

seeg = fiff.Raw('goumi-run3-seeg.fif', preload=True)


# Out[65]:

#     Opening raw data file /home/duke/tools/iidc/work/goumi-run3-seeg.fif...
#         Read a total of 1 projection items:
#             Average EEG reference (1 x 97)  idle
#         Range : 0 ... 1453885 =      0.000 ...   714.612 secs
#     Ready.
#     Reading 0 ... 1453885  =      0.000 ...   714.612 secs...
#     [done]
# 

# Reorganize for sEEG,

# In[87]:

import re

def seeg_ch_name_split(nm):
    elec, idx = re.match(r'([A-Za-z]+)(\d+)', nm).groups()
    return elec, int(idx)


# Sort with a custom key so 'A10' (=> 'A10') is after 'A2' (=> 'A02')

# In[67]:

chix = array([seeg.ch_names.index(nm) for nm in sorted(seeg.ch_names, key=lambda s: '%s%02d' % seeg_ch_name_split(nm))])
seeg.info['chs'] = [seeg.info['chs'][i] for i in chix]
seeg.info['ch_names'] = [seeg.info['ch_names'][i] for i in chix]
seeg._data = seeg._data[chix]


# In[70]:

seeg.save('goumi-run3-seeg.fif', overwrite=True, verbose=False)


# In[80]:

seeg = fiff.Raw('goumi-run3-seeg.fif', preload=True)


# Out[80]:

#     Opening raw data file /home/duke/tools/iidc/work/goumi-run3-seeg.fif...
#         Read a total of 1 projection items:
#             Average EEG reference (1 x 97)  idle
#         Range : 0 ... 1453885 =      0.000 ...   714.612 secs
#     Ready.
#     Reading 0 ... 1453885  =      0.000 ...   714.612 secs...
#     [done]
# 

# Check we've correctly sorted & don't have duplicate A channels:

# In[81]:

print seeg.ch_names[:12:3], [n for n in seeg.ch_names if n=='A1']


# Out[81]:

#     ['A1', 'A4', 'A7', 'A10'] ['A1']
# 

# I'm not sure how to go about setting up a bipolar montage. Maybe I should ask the question.
# 
# The direct way is just to rewrite channel names and info and form the data from the diff of `_data`, removing non-contiguous channels.

# In[88]:

def find_bip_idx(seeg):
    bip_idx = []
    for i in xrange(len(seeg.ch_names)-1):
        (e1, i1), (e2, i2) = map(seeg_ch_name_split, seeg.ch_names[i:i+2])
        if e1==e2:
            bip_idx.append((i, i+1, e1, i1, i2))
    return bip_idx


bip_data = []
bip_chan = []

#   16, 17, 'OF',   2,   3  results in 'OF2-3', _data[16] = _data[16] - _data[17]
for i1, i2, elec, ei1, ei2 in find_bip_idx(seeg):
    ch_info = seeg.info['chs'][i1].copy()
    ch_info['ch_name'] = '%s%d-%d' % (elec, ei1, ei2)
    bip_chan.append(ch_info)
    bip_data.append(seeg._data[i1] - seeg._data[i2])
    
seeg.info['chs'] = bip_chan
seeg._data = array(bip_data)
seeg.info['ch_names'] = [c['ch_name'] for c in seeg.info['chs']]
seeg.info['nchan'] = seeg._data.shape[0]

# check
seeg.ch_names[:12]


# Out[88]:

['A1-2',
 'A2-3',
 'A3-4',
 'A4-5',
 'A5-6',
 'A6-7',
 'A7-8',
 'A8-9',
 'A9-10',
 'A10-11',
 'B1-2',
 'B2-3']

util.fix_bti_import_channel_names is defined in a iidc.util module.