Filtering and Epoching MNE-HCP

External Email - Use Caution

Hi all,
I am still working with the MNE-HCP tools, and have two pipelines that I seem to have hit some walls for both of them. The overall goal is to commute coherence connectivity in source space.
In the first pipeline I am simply trying to use:
raw1=hcp.read_epochs(subject, 'task_working_memory', onset='stim', run_index=0, hcp_path=hcp_path, return_fixations_motor=False)
raw2=hcp.read_epochs(subject, 'task_working_memory', onset='stim', run_index=1, hcp_path=hcp_path, return_fixations_motor=False)
But I get:
line 611, in _read_epochs
    epochs.times = times
AttributeError: can't set attribute

Since I cant seem to get that working I have also been running the raw instance through my lab?s standard preprocessing pipeline, but there are some things in the tutorial that I am not sure if I have to include that are not part of our standard (eeg) pipeline. When I plot the sensors it appears 90 degrees off, but hcp.preprocessing.map_ch_coords_to_mne(raw) seems to fix that. However, I am not sure if I need to apply hcp.preprocessing.apply_ref_correction(raw) because I am getting some unexpected things by doing it. when I do:
ecg_evoked = mne.preprocessing.create_ecg_epochs(raw,ch_name='ECG-').average()
ecg_evoked.apply_baseline((None, None))
ecg_evoked.plot_joint()
projs, events = mne.preprocessing.compute_proj_ecg(raw, ch_name = 'ECG-', n_grad=1, n_mag=1, n_eeg=1, reject=None)
mne.viz.plot_projs_topomap(projs, info=raw.info)
raw.add_proj(projs)
without the reference correction it smears the cardiac artifact (not sure if I should be using ECG- or ECG+ or both) as a negative going component across the back left and front right of the head. With the reference corrections applied it gives me a very standard looking cardiac component but on the other side of the head near the left eye. I am use to EEG data but something odd I noticed is that there appears to be regular cardiac activity not just on the back left sensors, but also very strongly on both sides of the eyes, too regular for eyeblinks and matched up in time with the sensors that to typically pick up the cardiac artifact.

Does anyone out there know more of what is going on, I am at a loss and not sure what to do. Thank you for any advice!

Best,
Matt
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/attachments/20201001/4e2d1c78/attachment.html

External Email - Use Caution

Hi Matt,

Do you mind sharing a full script with me that I can use to debug?
It may turn out that we have a bit of maintenance to do in the mne-hcp
code to catch up with newer versions of main MNE etc.
Once we have found the problem I can fix it directly.

Best,
Denis

External Email - Use Caution

Of course, and thanks for looking into this. I have three pipelines the first one is short, just trying to read the epochs but gives me the epochs.times error:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import nibabel as nib
import pylab as pl
import mne
import scipy
from mne.time_frequency.tfr import cwt, morlet
import sys
sys.path.append('/N/u/mattwint/Carbonate/Desktop/net/')
import hcp
from hcp import preprocessing as preproc

subject='100307'
hcp_path='/N/dcwan/projects/hcp'
raw1=hcp.read_epochs(subject, 'task_working_memory', onset='stim', run_index=0, hcp_path=hcp_path, return_fixations_motor=False)
raw2=hcp.read_epochs(subject, 'task_working_memory', onset='stim', run_index=1, hcp_path=hcp_path, return_fixations_motor=False)

The second script is my attempt to work with the raw instance and take into account some of the things the GitHub mne-hcp page says to do, although I do not understand it all. I can run this script all they way through, but end up getting odd cardiac components I do not trust. I have tried running it with and without the reference correction and with the ecg+ and ecg- (but not both at the same time) and I have not gotten a component back that looks right to me, although I am more use to eeg data, so maybe it is just my ignorance. Note I pull the trial info, but have not gotten to the part were I epoch it yet:.
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import nibabel as nib
import pylab as pl
import mne
import scipy
from mne.time_frequency.tfr import cwt, morlet
import sys
sys.path.append('/N/u/mattwint/Carbonate/Desktop/data/')
import hcp

data_type='task_working_memory'
subject='100307'
hcp_path='/N/dcwan/projects/hcp'
raw1=hcp.read_raw(subject, data_type, run_index=0, hcp_path=hcp_path)
raw2=hcp.read_raw(subject, data_type, run_index=1, hcp_path=hcp_path)

#%%########################### trial info #####################################
trial_info1=hcp.read_trial_info(subject, data_type, run_index=0, hcp_path=hcp_path)
trial_info2=hcp.read_trial_info(subject, data_type, run_index=1, hcp_path=hcp_path)
info1=len(trial_info1['stim']['codes'])
info2=len(trial_info2['stim']['codes'])
trial_info=np.zeros([info1 + info2, 40])
trial_info[0:info1,0:40]=(trial_info1['stim']['codes'])
trial_info[info1:info1+info2,0:40]=(trial_info2['stim']['codes'])

#%%
concatenated=mne.concatenate_raws([raw1,raw2])
concatenated.annotations
concatenated.annotations.delete([0,1])
concatenated.plot()
raw = concatenated.copy()
raw.load_data()

#%%###################### ??? regress out ref. ??? ############################
hcp.preprocessing.apply_ref_correction(raw)

#%%#Filter the data to 0.1-100 Hz (IIR for HP, FIR for LP, as was established previously)
raw.filter(l_freq = 0.1, h_freq = None, method = 'iir')
raw.filter(h_freq = 100, l_freq = None, method = 'fir')

#%%#View the PSD to make sure our filter took effect and orient cap to head
hcp.preprocessing.map_ch_coords_to_mne(raw)
raw.plot_psd(area_mode='range', average=False, fmax = 200)

#%%
room_noise=hcp.read_raw(subject, 'noise_empty_room', hcp_path=hcp_path)
room_info = hcp.read_info(subject, 'noise_empty_room', hcp_path=hcp_path)

#%%
empty_room_projs = mne.compute_proj_raw(room_noise, n_grad=3, n_mag=3)
raw.add_proj(empty_room_projs)

#%%
system_projs = raw.info['projs']

mags = mne.pick_types(raw.info, meg='mag')
for title, projs in [('system', system_projs),
                     ('subject-specific', empty_room_projs[3:])]:
    raw.add_proj(projs, remove_existing=True)
    fig = raw.plot(proj=True, order=mags, duration=1, n_channels=2)
    fig.subplots_adjust(top=0.9) # make room for title
    fig.suptitle('{} projectors'.format(title), size='xx-large', weight='bold')

#%% ???(Do I need + & -)???
#The following few cells describe the process of deleting ECG projections
ecg_evoked = mne.preprocessing.create_ecg_epochs(raw,ch_name='ECG-').average()
ecg_evoked.apply_baseline((None, None))
ecg_evoked.plot_joint()

#Let's try making a projection based on our ECG data
projs, events = mne.preprocessing.compute_proj_ecg(raw, ch_name = 'ECG-', n_grad=1, n_mag=1, n_eeg=1, reject=None)

#Here we will print out the list of projectors
print(projs)

#This creates a scalp graphic with the distribution of projectors
mne.viz.plot_projs_topomap(projs, info=raw.info)
raw.add_proj(projs)

#%%#Let's do the same for vertical EOG
eog_evoked = mne.preprocessing.create_eog_epochs(raw, ch_name = 'VEOG-').average()
eog_evoked.apply_baseline((None, None))
eog_evoked.plot_joint()

#Compute projections; no_proj is set to "True" to avoid printing the whole set of projections computed above
eog_projs, _ = mne.preprocessing.compute_proj_eog(raw, n_grad=1, n_mag=1, n_eeg=1, reject=None, ch_name = 'VEOG-', no_proj=True)
#Plot the results of projections of vertical EOG
mne.viz.plot_projs_topomap(eog_projs, info=raw.info)
for title in ('Without', 'With'):
    if title == 'With':
        raw.add_proj(eog_projs)
    fig = raw.plot()
    fig.subplots_adjust(top=0.9) # make room for title
    fig.suptitle('{} EOG projectors'.format(title), size='xx-large',
                 weight='bold')

mne.viz.plot_projs_topomap(eog_projs, info=raw.info)
raw.add_proj(eog_projs)

#%%#And horizontal EOG
eogh_evoked = mne.preprocessing.create_eog_epochs(raw, ch_name = 'HEOG-').average()
eogh_evoked.apply_baseline((None, None))
eogh_evoked.plot_joint()

#Compute projections; no_proj is set to "True" to avoid printing the whole set of projections computed above
eogh_projs, _ = mne.preprocessing.compute_proj_eog(raw, n_grad=1, n_mag=1, n_eeg=1, reject=None, ch_name = 'HEOG-', no_proj=True)
#Plot the results of projections of horizontal EOG
mne.viz.plot_projs_topomap(eogh_projs, info=raw.info)
for title in ('Without', 'With'):
    if title == 'With':
        raw.add_proj(eogh_projs)
    fig = raw.plot()
    fig.subplots_adjust(top=0.9) # make room for title
    fig.suptitle('{} EOG projectors'.format(title), size='xx-large',
                weight='bold')

mne.viz.plot_projs_topomap(eogh_projs, info=raw.info)
raw.add_proj(eogh_projs)

#%%
raw.apply_proj

#%%#Let's back it up one more time
raw.save('concatenated'+'filt01-100_'+'SSPapplied_tsss.fif', overwrite=True)

#%%#Filter the data to 1-40, as was established previously
raw.filter(l_freq = 1, h_freq = None, method = 'iir')
raw.filter(h_freq = 55, l_freq = None, method = 'fir')

#%%#And back it up
raw.save('concatenated'+'filt1-40_tsss.fif', overwrite=True)

#%%
raw.plot()
raw.plot_psd(fmin = 1, fmax = 40)

My third script works fine, I pull over parts of the tmegpreproc files and concatenate them back together. The only problem I run into on that one is there are a different number of good sensors between run 1 and run 2, so I have yet to make a decision on the best way to work around that problem. If you want to see that I can send it too, it just doesn?t really utilize the mne-hcp tools.