Clean-up work

01_tsss.py

@@ -5,46 +5,29 @@
 
 @author: jiaxin
 """
-import os
-# os.chdir()
 import mne
 import argparse
-from mne_bids import BIDSPath, read_raw_bids, write_raw_bids, get_anonymization_daysback
+from mne_bids import BIDSPath, read_raw_bids, write_raw_bids
 from config import fname, bad_channels, event_id, task
 
 # %% Handle command line arguments
 parser = argparse.ArgumentParser(description=__doc__)
-parser.add_argument('--subject',  type=int, default=27,
-                    help='The subject to process', required=False)
-parser.add_argument('--runs',  type=int, default=1,
-                    help='how many runs', required=False)
-parser.add_argument('--run',  type=int, default=1,
-                    help='which run', required=False)
+parser.add_argument('--subject',  type=int, help='The subject to process')
+parser.add_argument('--run',  type=int, default=1, help='which run (subject 16 has 2 runs)', required=False)
 
 args = parser.parse_args()
 subject = args.subject
 print('Processing subject:', subject)
 
 bads = bad_channels[f'sub-{subject:02}']
-print('Bad channles:', bads)
+print('Bad channels:', bads)
 
 # mne.set_log_level(verbose="warning")
 tsss_log = fname.log(subject=subject, proc='tsss')
 #
 mne.set_log_file(tsss_log, overwrite=True)
 
-if args.runs == 1:
-    bids_path = BIDSPath(subject=f'{subject:02}', task=task,
-                         datatype='meg',
-                         root=fname.bids_dir
-                         )
-    write_bids_path = BIDSPath(subject=f'{subject:02}', task=task, datatype='meg',
-                               processing='tsss',
-                               root=fname.derivatives_dir)
-    raw = read_raw_bids(bids_path=bids_path, verbose=False)
-
-    dev_head_t_ref = raw.info['dev_head_t']
-else:
+if args.subject == 16:
     bids_path = BIDSPath(subject=f'{subject:02}', task=task,
                          datatype='meg',
                          run=f'{args.run:01}',
@@ -64,9 +47,18 @@
     raw = mne.io.read_raw_fif(bids_path, verbose=False)
     dev_head_t_ref = read_raw_bids(
         bids_path=bids_path1, verbose=False).info['dev_head_t']
+else:
+    bids_path = BIDSPath(subject=f'{subject:02}', task=task,
+                         datatype='meg',
+                         root=fname.bids_dir
+                         )
+    write_bids_path = BIDSPath(subject=f'{subject:02}', task=task, datatype='meg',
+                               processing='tsss',
+                               root=fname.derivatives_dir)
+    raw = read_raw_bids(bids_path=bids_path, verbose=False)
+
+    dev_head_t_ref = raw.info['dev_head_t']
 
-# raw.plot(lowpass=40)#%%
-# %%break annotationif
 if subject == 1:
     onsets = [raw.first_time+477, raw.first_time+1048]
     durations = [18, 54]
@@ -76,7 +68,6 @@
     )
     raw.set_annotations(raw.annotations + break_annot)
 
-#%%31.5.2024
 raw.info["bads"] = []
 auto_noisy_chs, auto_flat_chs, auto_scores = mne.preprocessing.find_bad_channels_maxwell(
          raw,
@@ -85,12 +76,11 @@
          return_scores=True,
          verbose=True,
      )
-raw.info["bads"] = list(set(bads+ auto_noisy_chs + auto_flat_chs))
-    
-print(raw.info["bads"])
+raw.info["bads"] = list(set(bads + auto_noisy_chs + auto_flat_chs))
+
+print("Bad channels after maxfilter detection:", raw.info["bads"])
 
 
-# %%
 chpi_amplitudes = mne.chpi.compute_chpi_amplitudes(raw, verbose=True)
 chpi_locs = mne.chpi.compute_chpi_locs(raw.info, chpi_amplitudes, verbose=True)
 head_pos = mne.chpi.compute_head_pos(raw.info, chpi_locs, verbose=True)
@@ -102,15 +92,13 @@
     st_duration=60,
     coord_frame='head',
     destination=dev_head_t_ref,
-
     head_pos=head_pos,
     verbose=True
 )
 raw_tsss = mne.chpi.filter_chpi(raw_tsss)
 raw_tsss.drop_channels([f'CHPI{i:03d}' for i in range(1, 10)])
 
 
-# %%%
 event_id["BAD_ACQ_SKIP"] = 0
 write_raw_bids(
     raw=raw_tsss,

02_filter.py

@@ -7,15 +7,13 @@
 import mne
 import argparse
 from mne_bids import BIDSPath, read_raw_bids, write_raw_bids
-from config import fname, bandpass_fmin, bandpass_fmax, task, bad_channels, phase, event_id
+from config import fname, bandpass_fmin, bandpass_fmax, task, event_id
 
 
 # %%
 parser = argparse.ArgumentParser(description=__doc__)
-parser.add_argument('--subject',  type=int, default=27,
-                    help='The subject to process')
-parser.add_argument('--run',  type=int, default=None,
-                    help='which run', required=False)
+parser.add_argument('--subject',  type=int, help='The subject to process')
+parser.add_argument('--run',  type=int, default=1, help='which run (subject 16 has 2 runs)', required=False)
 args = parser.parse_args()
 subject = args.subject
 print('Processing subject:', subject)

03_ica.py

@@ -8,12 +8,12 @@
 import argparse
 import numpy as np
 import mne
-from mne.preprocessing import ICA, create_ecg_epochs, create_eog_epochs, find_ecg_events, find_eog_events
+from mne.preprocessing import ICA, create_eog_epochs, find_ecg_events
 from mne_bids import BIDSPath, read_raw_bids
-from config import (fname, n_ecg_components,
-                    n_eog_components, bad_channels, task, eog_chs)
-from mne.io.pick import _picks_to_idx, pick_types, pick_channels
+from config import (fname, n_ecg_components, n_eog_components, task)
+from mne.io.pick import _picks_to_idx
 from mne.epochs import Epochs
+
 # %%
 # Be verbose
 mne.set_log_level('INFO')
@@ -26,6 +26,7 @@
 args = parser.parse_args()
 subject = f'sub-{args.subject:02}'
 print('Processing subject:', subject)
+
 # %%
 run = None if args.run == None else f'{args.run:01}'
 # Construct a raw object that will load the highpass-filtered data.
@@ -36,9 +37,6 @@
 raw = read_raw_bids(bids_path=bids_path, verbose=False)
 
 # %%
-# raw.crop(tmax=578) #subject 12 EOG1 didn't work from 578 s
-# raw.crop(tmax=455) #subject 1
-# %% fitting ica
 n_components = 0.99
 print('Fitting ICA')
 ica = ICA(method='fastica', random_state=42, n_components=n_components)
@@ -74,16 +72,7 @@
     0.1, 40, l_trans_bandwidth='auto',
     h_trans_bandwidth='auto', filter_length='auto', phase='zero',
     fir_window='hamming', fir_design='firwin', n_jobs=-1)
-# %% find ecg events for plotting
-# import matplotlib.pyplot as plt
-# events, a, b, ecg = find_ecg_events(
-#         raw2,
-#         l_freq=8,
-#         h_freq=16,
-#         return_ecg=True,
-#         reject_by_annotation=True,
-#     )
-# %%
+
 event_id = 999
 ecg_events, _, _, ecg = find_ecg_events(
     raw2,

04_epochs.py

@@ -17,18 +17,21 @@
 mne.set_log_level('INFO')
 
 parser = argparse.ArgumentParser(description=__doc__)
-parser.add_argument('--subject',  type=int, default=16,
-                    help='The subject to process', required=False)
-parser.add_argument('--runs',  type=int, default=1,
-                    help='how many runs', required=False)
+parser.add_argument('--subject',  type=int, help='The subject to process')
 args = parser.parse_args()
 subject = f'sub-{args.subject:02}'
 print('Processing subject:', subject)
 print(fname.report_html(subject=subject))
 # %%combine epochs
 epochs_list = []
-for run in range(args.runs):
-    run = f'{run+1:01}' if args.runs > 1 else None
+
+if subject == "sub-16":
+    n_runs = 2
+else:
+    n_runs = 1
+
+for run in range(n_runs):
+    run = f'{run+1:01}' if n_runs > 1 else None
 # Construct a raw object that will load the highpass-filtered data.
     bids_path = BIDSPath(subject=subject[-2:], task=task, datatype='meg', processing='filt',
                          run=run,

05_visualize.py

@@ -27,8 +27,7 @@
 mne.set_log_level('INFO')
 
 parser = argparse.ArgumentParser(description=__doc__)
-parser.add_argument('--subject',  type=int, default=1,
-                    help='The subject to process', required=False)
+parser.add_argument('--subject',  type=int, help='The subject to process')
 args = parser.parse_args()
 subject = f'sub-{args.subject:02}'
 print('Processing subject:', subject)

07_subject_inv.py

@@ -6,20 +6,14 @@
 # %% import modules
 import argparse
 import mne
-from mne.minimum_norm import make_inverse_operator, apply_inverse, write_inverse_operator, read_inverse_operator
-from config import fname, event_id, spacing
-import os
-import os.path as op
-import matplotlib.pyplot as plt
-from mne.datasets import fetch_fsaverage
-from mpl_toolkits.axes_grid1 import make_axes_locatable
+from mne.minimum_norm import make_inverse_operator, write_inverse_operator
+from config import fname, spacing
 # %%
 # Be verbose
 mne.set_log_level('INFO')
 
 parser = argparse.ArgumentParser(description=__doc__)
-parser.add_argument('--subject',  type=int, default=1,
-                    help='The subject to process', required=False)
+parser.add_argument('--subject',  type=int, help='The subject to process')
 args = parser.parse_args()
 
 # MEG-MRI co-registration??
@@ -31,90 +25,37 @@
 
 # %% Read the epochs
 print('Reading epochs...')
-epochs = mne.read_epochs(
-    fname.epo(subject=subject))
+epochs = mne.read_epochs(fname.epo(subject=subject))
 
-# %% Prepare brain
-# epochs.pick_types('grad')
-# ave = {cat: epochs[cat].average() for cat in event_id}
-# %% config
-SUBJECT = subject
 mne.set_config('SUBJECTS_DIR', fname.mri_subjects_dir)
-# %%
-# MEG-MRI co-registration
-# if co_reg:
-# mne.gui.coregistration(subject=SUBJECT, inst=fname.epo(
-#     subject=subject),
-#     # trans=transname
-# )
-# %%
-# empty-room covariance, which captures noise from the sensors and environment.
+
 noise_cov = mne.compute_covariance(
     epochs, tmax=0.0, method='auto', rank='info')
 noise_cov = mne.cov.regularize(noise_cov, epochs.info)
-# noise_cov.plot(epochs.info)
-# %%
-# noise_cov = mne.compute_covariance(epochs, method='shrunk', rank='info')
-# noise_cov = mne.cov.regularize(noise_cov, epochs.info)
-# cov.plot()
-
-# %%
-# # Create the BEM model using the FreeSurfer watershed algorithm
-# bem=mne.bem.make_watershed_bem(
-#     subject=SUBJECT, atlas=True
-#     )
-
 
 # %% Establish the source space
-
-src = mne.setup_source_space(subject=SUBJECT, spacing=spacing, add_dist=False)
-src.save(fname.src(subject=SUBJECT), overwrite=True)
+src = mne.setup_source_space(subject=subject, spacing=spacing, add_dist=False)
+src.save(fname.src(subject=subject), overwrite=True)
 
 
 # %% Create the BEM model
 model_surfaces = mne.make_bem_model(
-    subject=SUBJECT, conductivity=[0.3])
+    subject=subject, conductivity=[0.3])
 bem = mne.make_bem_solution(model_surfaces)
 
 
 # Compute the lead-field matrix, i.e, the forward solution
-
-
-# if os.path.isfile(fwdname):
-#     fwd = mne.read_forward_solution(fwdname)
-# else:
 fwd = mne.make_forward_solution(
-    epochs.info, trans=fname.trans(subject=SUBJECT), src=src, bem=bem, n_jobs=-1,
+    epochs.info, trans=fname.trans(subject=subject), src=src, bem=bem, n_jobs=-1,
     # mindist=5,
     eeg=False)
-mne.write_forward_solution(fname.fwd_r(subject=SUBJECT), fwd, overwrite=True)
+mne.write_forward_solution(fname.fwd_r(subject=subject), fwd, overwrite=True)
 fwd = mne.convert_forward_solution(fwd, surf_ori=True)
 # Use only MEG channels
 fwd = mne.pick_types_forward(fwd, meg=True, eeg=False)
 
-# # To avoid re-computing fwd every time, you can d
-
-# %%
 inv = make_inverse_operator(epochs.info, fwd, noise_cov,
                             loose=0.2,
                             # loose='auto',
                             depth=0.8)
-write_inverse_operator(fname.inv(subject=SUBJECT), inv, overwrite=True)
-
-# %%
-# mne.viz.plot_alignment(
-#     info=epochs.info,
-#     trans=fname.trans(subject=SUBJECT),
-#     subject=SUBJECT,
-#     src=src,
-#     # subjects_dir=subjects_dir,
-#     # dig=True,
-#     # surfaces=["head-dense", "white"],
-#     # coord_frame="meg",
-# )
-# surfaces = dict(brain=0.4, outer_skull=0.6, head=None)
-# mne.viz.plot_alignment(info=epochs.info, trans=transname,
-#                        subject=SUBJECT, surfaces=surfaces,
-#                        fwd=fwd, src=src,
-#                        mri_fiducials=True,
-#                        bem=bem)
+write_inverse_operator(fname.inv(subject=subject), inv, overwrite=True)