diff --git a/examples/camcan/preprocess.py b/examples/camcan/1_preprocess.py similarity index 72% rename from examples/camcan/preprocess.py rename to examples/camcan/1_preprocess.py index d34e7784..7e297cde 100644 --- a/examples/camcan/preprocess.py +++ b/examples/camcan/1_preprocess.py @@ -2,27 +2,33 @@ """ +# Authors: Chetan Gohil + import pathlib from glob import glob from dask.distributed import Client + from osl import preprocessing, utils +# Directories BASE_DIR = "/well/woolrich/projects/camcan" RAW_DIR = BASE_DIR + "/cc700/meg/pipeline/release005/BIDSsep/derivatives_rest/aa/AA_movecomp/aamod_meg_maxfilt_00002" -PREPROC_DIR = BASE_DIR + "/winter23/preproc" +PREPROC_DIR = BASE_DIR + "/summer23/preproc" +# Settings config = """ preproc: - crop: {tmin: 30} - - filter: {l_freq: 0.5, h_freq: 125, method: iir, iir_params: {order: 5, ftype: butter}} - - notch_filter: {freqs: 50 100} - - resample: {sfreq: 250} + - filter: {l_freq: 0.5, h_freq: 200, method: iir, iir_params: {order: 5, ftype: butter}} + - notch_filter: {freqs: 50 88 100 118 150 156 185, notch_widths: 2} + - resample: {sfreq: 400} - bad_segments: {segment_len: 500, picks: mag, significance_level: 0.1} - bad_segments: {segment_len: 500, picks: grad, significance_level: 0.1} - bad_segments: {segment_len: 500, picks: mag, mode: diff, significance_level: 0.1} - bad_segments: {segment_len: 500, picks: grad, mode: diff, significance_level: 0.1} - - bad_channels: {picks: meg, significance_level: 0.1} - - ica_raw: {picks: meg, n_components: 64} + - bad_channels: {picks: mag, significance_level: 0.1} + - bad_channels: {picks: grad, significance_level: 0.1} + - ica_raw: {picks: meg, n_components: 0.99} - ica_autoreject: {picks: meg, ecgmethod: correlation, eogthreshold: auto} - interpolate_bads: {} """ @@ -32,7 +38,7 @@ # Get input files inputs = [] - for subject in glob(RAW_DIR + "/sub-*"): + for subject in sorted(glob(RAW_DIR + "/sub-*")): subject = pathlib.Path(subject).stem inputs.append(RAW_DIR + f"/{subject}/mf2pt2_{subject}_ses-rest_task-rest_meg.fif") diff --git a/examples/camcan/2_coregister.py b/examples/camcan/2_coregister.py new file mode 100644 index 00000000..9702f731 --- /dev/null +++ b/examples/camcan/2_coregister.py @@ -0,0 +1,108 @@ +"""Coregisteration. + +The scripts was first run for all subjects (with n_init=1). Then for subjects +whose coregistration looked a bit off we re-run this script just for that +particular subject with a higher n_init. +""" + +# Authors: Chetan Gohil + +import numpy as np +import pathlib +from glob import glob +from dask.distributed import Client + +from osl import source_recon, utils + +# Directories +BASE_DIR = "/well/woolrich/projects/camcan" +PREPROC_DIR = BASE_DIR + "/summer23/preproc" +COREG_DIR = BASE_DIR + "/summer23/coreg" +ANAT_DIR = BASE_DIR + "/cc700/mri/pipeline/release004/BIDS_20190411/anat" +FSL_DIR = "/well/woolrich/projects/software/fsl" + +# Files +PREPROC_FILE = ( + PREPROC_DIR + + "/mf2pt2_{0}_ses-rest_task-rest_meg" + + "/mf2pt2_{0}_ses-rest_task-rest_meg_preproc_raw.fif" +) +SMRI_FILE = ANAT_DIR + "/{0}/anat/{0}_T1w.nii" + +# Settings +config = """ + source_recon: + - extract_fiducials_from_fif: {} + - fix_headshape_points: {} + - compute_surfaces: + include_nose: False + - coregister: + use_nose: False + use_headshape: True + #n_init: 50 +""" + + +def fix_headshape_points(src_dir, subject, preproc_file, smri_file, epoch_file): + filenames = source_recon.rhino.get_coreg_filenames(src_dir, subject) + + # Load saved headshape and nasion files + hs = np.loadtxt(filenames["polhemus_headshape_file"]) + nas = np.loadtxt(filenames["polhemus_nasion_file"]) + lpa = np.loadtxt(filenames["polhemus_lpa_file"]) + rpa = np.loadtxt(filenames["polhemus_rpa_file"]) + + # Drop nasion by 4cm and remove headshape points more than 7 cm away + nas[2] -= 40 + distances = np.sqrt( + (nas[0] - hs[0]) ** 2 + (nas[1] - hs[1]) ** 2 + (nas[2] - hs[2]) ** 2 + ) + keep = distances > 70 + hs = hs[:, keep] + + # Remove anything outside of rpa + keep = hs[0] < rpa[0] + hs = hs[:, keep] + + # Remove anything outside of lpa + keep = hs[0] > lpa[0] + hs = hs[:, keep] + + # Remove headshape points on the neck + remove = hs[2] < min(lpa[2], rpa[2]) - 4 + hs = hs[:, ~remove] + + # Overwrite headshape file + utils.logger.log_or_print(f"overwritting {filenames['polhemus_headshape_file']}") + np.savetxt(filenames["polhemus_headshape_file"], hs) + + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + source_recon.setup_fsl(FSL_DIR) + + # Get subjects + subjects = [] + for subject in sorted(glob(PREPROC_FILE.replace("{0}", "*"))): + subjects.append(pathlib.Path(subject).stem.split("_")[1]) + + # Setup files + smri_files = [] + preproc_files = [] + for subject in subjects: + smri_files.append(SMRI_FILE.format(subject)) + preproc_files.append(PREPROC_FILE.format(subject)) + + # Setup parallel processing + # client = Client(n_workers=8, threads_per_worker=1) + + # Run coregistration + source_recon.run_src_batch( + config, + src_dir=COREG_DIR, + subjects=subjects, + preproc_files=preproc_files, + smri_files=smri_files, + extra_funcs=[fix_headshape_points], + # dask_client=True, + ) diff --git a/examples/camcan/source_reconstruct.py b/examples/camcan/3_source_reconstruct.py similarity index 63% rename from examples/camcan/source_reconstruct.py rename to examples/camcan/3_source_reconstruct.py index e34ca77b..8aac6ffd 100644 --- a/examples/camcan/source_reconstruct.py +++ b/examples/camcan/3_source_reconstruct.py @@ -7,27 +7,35 @@ import pathlib from glob import glob from dask.distributed import Client + from osl import source_recon, utils +# Directories BASE_DIR = "/well/woolrich/projects/camcan" -PREPROC_DIR = BASE_DIR + "/winter23/preproc" -SRC_DIR = BASE_DIR + "/winter23/src" -ANAT_DIR = BASE_DIR + "/cc700/mri/pipeline/release004/BIDS_20190411/anat" -PREPROC_FILE = PREPROC_DIR + "/mf2pt2_{0}_ses-rest_task-rest_meg/mf2pt2_{0}_ses-rest_task-rest_meg_preproc_raw.fif" -SMRI_FILE = ANAT_DIR + "/{0}/anat/{0}_T1w.nii" +PREPROC_DIR = BASE_DIR + "/summer23/preproc" +SRC_DIR = BASE_DIR + "/summer23/src" FSL_DIR = "/well/woolrich/projects/software/fsl" +# Files +PREPROC_FILE = ( + PREPROC_DIR + + "/mf2pt2_{0}_ses-rest_task-rest_meg" + + "/mf2pt2_{0}_ses-rest_task-rest_meg_preproc_raw.fif" +) + +# Settings config = """ source_recon: - forward_model: model: Single Layer - beamform_and_parcellate: - freq_range: [1, 45] + freq_range: [1, 80] chantypes: [mag, grad] rank: {meg: 60} - parcellation_file: fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz + parcellation_file: Glasser52_binary_space-MNI152NLin6_res-8x8x8.nii.gz method: spatial_basis orthogonalisation: symmetric + extra_chans: [eog, ecg] """ if __name__ == "__main__": @@ -36,20 +44,12 @@ # Get subjects subjects = [] - for subject in sorted( - glob( - PREPROC_DIR - + "/mf2pt2_*_ses-rest_task-rest_meg" - + "/mf2pt2_sub-*_ses-rest_task-rest_meg_preproc_raw.fif" - ) - ): + for subject in sorted(glob(PREPROC_FILE.replace("{0}", "*"))): subjects.append(pathlib.Path(subject).stem.split("_")[1]) # Setup files - smri_files = [] preproc_files = [] for subject in subjects: - smri_files.append(SMRI_FILE.format(subject)) preproc_files.append(PREPROC_FILE.format(subject)) # Setup parallel processing @@ -61,6 +61,5 @@ src_dir=SRC_DIR, subjects=subjects, preproc_files=preproc_files, - smri_files=smri_files, dask_client=True, ) diff --git a/examples/camcan/sign_flip.py b/examples/camcan/4_sign_flip.py similarity index 86% rename from examples/camcan/sign_flip.py rename to examples/camcan/4_sign_flip.py index 210d0f0c..c3b805f6 100644 --- a/examples/camcan/sign_flip.py +++ b/examples/camcan/4_sign_flip.py @@ -2,12 +2,16 @@ """ +# Authors: Chetan Gohil + from glob import glob from dask.distributed import Client + from osl import utils from osl.source_recon import find_template_subject, run_src_batch, setup_fsl -SRC_DIR = "/well/woolrich/projects/camcan/winter23/src" +# Directories +SRC_DIR = "/well/woolrich/projects/camcan/summer23/src" FSL_DIR = "/well/woolrich/projects/software/fsl" if __name__ == "__main__": @@ -30,8 +34,8 @@ template: {template} n_embeddings: 15 standardize: True - n_init: 3 - n_iter: 3000 + n_init: 5 + n_iter: 5000 max_flips: 20 """ diff --git a/examples/camcan/5_prepare.py b/examples/camcan/5_prepare.py new file mode 100644 index 00000000..8ecb004b --- /dev/null +++ b/examples/camcan/5_prepare.py @@ -0,0 +1,19 @@ +"""Data preparation: time-delay embedding and principal component analysis. + +""" + +# Authors: Chetan Gohil + +from glob import glob + +from osl_dynamics.data import Data + +files = sorted(glob("/well/woolrich/projects/camcan/summer23/src/*/sflip_parc-raw.fif")) +data = Data(files, picks="misc", reject_by_annotation="omit", n_jobs=16) +methods = { + "tde_pca": {"n_embeddings": 15, "n_pca_components": 100}, + "standardize": {}, +} +data.prepare(methods) +data.save("/well/woolrich/projects/camcan/summer23/prepared") +data.delete_dir() diff --git a/examples/camcan/calc_glm.py b/examples/camcan/calc_glm.py deleted file mode 100644 index f9bfd96b..00000000 --- a/examples/camcan/calc_glm.py +++ /dev/null @@ -1,196 +0,0 @@ -"""Calculate first-level GLM. - -""" - -import os -import osl -import mne -import h5py -import sails -import numpy as np -from glob import glob -from scipy import stats -from dask.distributed import Client - - -def get_device_fids(raw): - # Put fiducials in device space - head_fids = mne.viz._3d._fiducial_coords(raw.info["dig"]) - head_fids = np.vstack(([0, 0, 0], head_fids)) - fid_space = raw.info["dig"][0]["coord_frame"] - assert fid_space == 4 # Ensure we have FIFFV_COORD_HEAD coords - - # Get device to head transform and inverse - dev2head = raw.info["dev_head_t"] - head2dev = mne.transforms.invert_transform(dev2head) - assert head2dev["from"] == 4 - assert head2dev["to"] == 1 - - # Apply transformation to get fids in device space - device_fids = mne.transforms.apply_trans(head2dev, head_fids) - - return device_fids - - -def make_bads_regressor(raw, mode="raw"): - bads = np.zeros((raw.n_times,)) - for an in raw.annotations: - if an["description"].startswith("bad") and an["description"].endswith(mode): - start = raw.time_as_index(an["onset"])[0] - raw.first_samp - duration = int(an["duration"] * raw.info["sfreq"]) - bads[start : start + duration] = 1 - if mode == "raw": - bads[: int(raw.info["sfreq"] * 2)] = 1 - bads[-int(raw.info["sfreq"] * 2) :] = 1 - else: - bads[: int(raw.info["sfreq"])] = 1 - bads[-int(raw.info["sfreq"]) :] = 1 - return bads - - -def run_first_level(dataset, userargs): - run_id = osl.utils.find_run_id(dataset["raw"].filenames[0]) - subj_id = run_id.split("_")[0] - - # Bad segments - bads = make_bads_regressor(dataset["raw"], mode="mag") - - # EOGs - vertical only - eogs = dataset["raw"].copy().pick_types(meg=False, eog=True) - eogs = eogs.filter(l_freq=1, h_freq=20, picks="eog").get_data() - - # ECG - lots of sanity checking - ecg_events, ch_ecg, av_pulse = mne.preprocessing.find_ecg_events(dataset["raw"]) - ecg_events[:, 0] = ecg_events[:, 0] - dataset["raw"].first_samp - ecg = np.zeros_like(bads) - median_beat = np.median(np.diff(ecg_events[:, 0])) - last_beat = median_beat - for ii in range(ecg_events.shape[0] - 1): - beat = ecg_events[ii + 1, 0] - ecg_events[ii, 0] - if np.abs(beat - last_beat) > 50: - beat = last_beat - ecg[ecg_events[ii, 0] : ecg_events[ii + 1, 0]] = beat - beat = last_beat - ecg[ecg == 0] = median_beat - ecg = ecg / dataset["raw"].info["sfreq"] * 60 - - # Store covariates - confs = {"VEOG": np.abs(eogs[0, :]), "BadSegs": bads} - covs = {"Linear": np.linspace(-1, 1, ecg.shape[0]), "ECG": ecg} - conds = None - conts = None - - #%% ------------------------------------------- - # GLM - With z-transform - - # Get data - XX = stats.zscore(dataset["raw"].get_data(picks="meg")) - fs = dataset["raw"].info["sfreq"] - - # Fit model - freq_vect, copes, varcopes, extras = sails.stft.glm_periodogram( - XX, - axis=1, - fit_constant=True, - conditions=conds, - covariates=covs, - confounds=confs, - contrasts=conts, - nperseg=int(fs * 2), - noverlap=int(fs), - fmin=0.1, - fmax=100, - fs=fs, - mode="magnitude", - fit_method="glmtools", - ) - model, design, data = extras - - # Save out - outdir = userargs.get("outdir") - hdfname = os.path.join( - outdir, "{subj_id}-glm-data-ztrans.hdf5".format(subj_id=run_id) - ) - if os.path.exists(hdfname): - print("Overwriting previous results") - os.remove(hdfname) - with h5py.File(hdfname, "w") as F: - model.to_hdf5(F.create_group("model")) - design.to_hdf5(F.create_group("design")) - # data.to_hdf5(F.create_group('data')) - F.create_dataset("freq_vect", data=freq_vect) - F.create_dataset("device_fids", data=get_device_fids(dataset["raw"])) - F.create_dataset("scan_duration", data=dataset["raw"].times[-1]) - F.create_dataset("av_bpm", data=av_pulse) - - #%% ------------------------------------------- - # GLM - No z-transform - - # Get data - XX = dataset["raw"].get_data(picks="meg") - fs = dataset["raw"].info["sfreq"] - - # Fit model - freq_vect, copes, varcopes, extras = sails.stft.glm_periodogram( - XX, - axis=1, - fit_constant=True, - conditions=conds, - covariates=covs, - confounds=confs, - contrasts=conts, - nperseg=int(fs * 2), - noverlap=int(fs), - fmin=0.1, - fmax=100, - fs=fs, - mode="magnitude", - fit_method="glmtools", - ) - model, design, data = extras - - # Save out - outdir = userargs.get("outdir") - hdfname = os.path.join( - outdir, "{subj_id}-glm-data-notrans.hdf5".format(subj_id=run_id) - ) - if os.path.exists(hdfname): - print("Overwriting previous results") - os.remove(hdfname) - with h5py.File(hdfname, "w") as F: - model.to_hdf5(F.create_group("model")) - design.to_hdf5(F.create_group("design")) - # data.to_hdf5(F.create_group('data')) - F.create_dataset("freq_vect", data=freq_vect) - F.create_dataset("device_fids", data=get_device_fids(dataset["raw"])) - F.create_dataset("scan_duration", data=dataset["raw"].times[-1]) - F.create_dataset("av_bpm", data=av_pulse) - - # Always need to return the dataset - even though its unchanged in this func - return dataset - - -config = """ - preproc: - - run_first_level: {outdir: /well/woolrich/projects/camcan/winter23/glm} -""" - -if __name__ == "__main__": - inputs = sorted( - glob( - "/well/woolrich/projects/camcan/winter23/preproc" - + "/mf2pt2_*_ses-rest_task-rest_meg" - + "/mf2pt2_*_ses-rest_task-rest_meg_preproc_raw.fif" - ) - ) - - client = Client(n_workers=16, threads_per_worker=1) - - osl.preprocessing.run_proc_batch( - config, - inputs, - outdir="./tmp", # this directory can be deleted after running - extra_funcs=[run_first_level], - overwrite=True, - dask_client=True, - ) diff --git a/examples/camcan/coregister.py b/examples/camcan/coregister.py deleted file mode 100644 index f18fbeb9..00000000 --- a/examples/camcan/coregister.py +++ /dev/null @@ -1,273 +0,0 @@ -"""Coregisteration. - -""" - -import numpy as np -import pathlib -from glob import glob -from dask.distributed import Client -from osl import source_recon, utils - -BASE_DIR = "/well/woolrich/projects/camcan" -PREPROC_DIR = BASE_DIR + "/winter23/preproc" -COREG_DIR = BASE_DIR + "/winter23/coreg" -ANAT_DIR = BASE_DIR + "/cc700/mri/pipeline/release004/BIDS_20190411/anat" -PREPROC_FILE = PREPROC_DIR + "/mf2pt2_{0}_ses-rest_task-rest_meg/mf2pt2_{0}_ses-rest_task-rest_meg_preproc_raw.fif" -SMRI_FILE = ANAT_DIR + "/{0}/anat/{0}_T1w.nii" -FSL_DIR = "/well/woolrich/projects/software/fsl" - -config = """ - source_recon: - - extract_fiducials_from_fif: {} - - fix_headshape_points: {} - - custom_coregister: {} -""" - -def fix_headshape_points(src_dir, subject, preproc_file, smri_file, epoch_file): - filenames = source_recon.rhino.get_coreg_filenames(src_dir, subject) - - # Load saved headshape and nasion files - hs = np.loadtxt(filenames["polhemus_headshape_file"]) - nas = np.loadtxt(filenames["polhemus_nasion_file"]) - lpa = np.loadtxt(filenames["polhemus_lpa_file"]) - rpa = np.loadtxt(filenames["polhemus_rpa_file"]) - - # Drop nasion by 4cm - nas[2] -= 40 - distances = np.sqrt( - (nas[0] - hs[0]) ** 2 + (nas[1] - hs[1]) ** 2 + (nas[2] - hs[2]) ** 2 - ) - - # Keep headshape points more than 7 cm away - keep = distances > 70 - hs = hs[:, keep] - - # Remove anything outside of rpa - keep = hs[0] < rpa[0] - hs = hs[:, keep] - - # Remove anything outside of lpa - keep = hs[0] > lpa[0] - hs = hs[:, keep] - - if subject in [ - "sub-CC110606", "sub-CC120061", "sub-CC120727", "sub-CC221648", "sub-CC221775", - "sub-CC121397", "sub-CC121795", "sub-CC210172", "sub-CC210519", "sub-CC210657", - "sub-CC220419", "sub-CC220974", "sub-CC221107", "sub-CC222956", "sub-CC221828", - "sub-CC222185", "sub-CC222264", "sub-CC222496", "sub-CC310224", "sub-CC310361", - "sub-CC310414", "sub-CC312222", "sub-CC320022", "sub-CC320088", "sub-CC320321", - "sub-CC320336", "sub-CC320342", "sub-CC320448", "sub-CC322186", "sub-CC410243", - "sub-CC420091", "sub-CC420094", "sub-CC420143", "sub-CC420167", "sub-CC420241", - "sub-CC420261", "sub-CC420383", "sub-CC420396", "sub-CC420435", "sub-CC420493", - "sub-CC420566", "sub-CC420582", "sub-CC420720", "sub-CC510255", "sub-CC510321", - "sub-CC510323", "sub-CC510415", "sub-CC510480", "sub-CC520002", "sub-CC520011", - "sub-CC520042", "sub-CC520055", "sub-CC520078", "sub-CC520127", "sub-CC520239", - "sub-CC520254", "sub-CC520279", "sub-CC520377", "sub-CC520391", "sub-CC520477", - "sub-CC520480", "sub-CC520552", "sub-CC520775", "sub-CC610050", "sub-CC610576", - "sub-CC620354", "sub-CC620406", "sub-CC620479", "sub-CC620518", "sub-CC620557", - "sub-CC620572", "sub-CC620610", "sub-CC620659", "sub-CC621011", "sub-CC621128", - "sub-CC621642", "sub-CC710131", "sub-CC710350", "sub-CC710551", "sub-CC710591", - "sub-CC710982", "sub-CC711128", "sub-CC720119", "sub-CC720188", "sub-CC720238", - "sub-CC720304", "sub-CC720358", "sub-CC720511", "sub-CC720622", "sub-CC720685", - "sub-CC721292", "sub-CC721504", "sub-CC721519", "sub-CC721891", "sub-CC721374", - "sub-CC722542", "sub-CC722891", "sub-CC121111", "sub-CC121144", "sub-CC210250", - "sub-CC210422", "sub-CC220519", "sub-CC221209", "sub-CC221487", "sub-CC221595", - "sub-CC221886", "sub-CC310331", "sub-CC410121", "sub-CC410179", "sub-CC420157", - "sub-CC510395", "sub-CC610653", - ]: - # Remove headshape points 1cm below lpa - keep = hs[2] > (lpa[2] - 10) - hs = hs[:, keep] - elif subject in [ - "sub-CC210023", "sub-CC210124", "sub-CC220132", "sub-CC220203", "sub-CC220223", - "sub-CC221054", "sub-CC310410", "sub-CC320089", "sub-CC320651", "sub-CC721114", - "sub-CC320680", "sub-CC320759", "sub-CC320776", "sub-CC320850", "sub-CC320888", - "sub-CC320904", "sub-CC321073", "sub-CC321154", "sub-CC321174", "sub-CC410084", - "sub-CC410101", "sub-CC410173", "sub-CC410226", "sub-CC410287", "sub-CC410325", - "sub-CC410432", "sub-CC420089", "sub-CC420149", "sub-CC420197", "sub-CC420198", - "sub-CC420217", "sub-CC420222", "sub-CC420260", "sub-CC420324", "sub-CC420356", - "sub-CC420454", "sub-CC420589", "sub-CC420888", "sub-CC510039", "sub-CC510115", - "sub-CC510161", "sub-CC510237", "sub-CC510258", "sub-CC510355", "sub-CC510438", - "sub-CC510551", "sub-CC510609", "sub-CC510629", "sub-CC510648", "sub-CC512003", - "sub-CC520013", "sub-CC520065", "sub-CC520083", "sub-CC520097", "sub-CC520147", - "sub-CC520168", "sub-CC520209", "sub-CC520211", "sub-CC520215", "sub-CC520247", - "sub-CC520395", "sub-CC520398", "sub-CC520503", "sub-CC520584", "sub-CC610052", - "sub-CC610178", "sub-CC610210", "sub-CC610212", "sub-CC610288", "sub-CC610575", - "sub-CC610631", "sub-CC620085", "sub-CC620090", "sub-CC620121", "sub-CC620164", - "sub-CC620444", "sub-CC620496", "sub-CC620526", "sub-CC620592", "sub-CC620793", - "sub-CC620919", "sub-CC710313", "sub-CC710486", "sub-CC710566", "sub-CC720023", - "sub-CC720497", "sub-CC720516", "sub-CC720646", "sub-CC721224", "sub-CC721729", - "sub-CC723395", "sub-CC222326", "sub-CC310160", "sub-CC121479", "sub-CC121685", - "sub-CC221755", "sub-CC320687", "sub-CC620152", "sub-CC711244", - ]: - # Remove headshape points below rpa - keep = hs[2] > rpa[2] - hs = hs[:, keep] - elif subject in [ - "sub-CC210617", "sub-CC220107", "sub-CC220198", "sub-CC220234", "sub-CC220323", - "sub-CC220335", "sub-CC222125", "sub-CC222258", "sub-CC310008", "sub-CC610046", - "sub-CC610508", - ]: - # Remove headshape points on the face - keep = np.logical_or(hs[2] > lpa[2], hs[1] < lpa[1]) - hs = hs[:, keep] - elif subject in [ - "sub-CC410129", "sub-CC410222", "sub-CC410323", "sub-CC410354", "sub-CC420004", - "sub-CC410390", "sub-CC420348", "sub-CC420623", "sub-CC420729", "sub-CC510043", - "sub-CC510086", "sub-CC510304", "sub-CC510474", "sub-CC520122", "sub-CC521040", - "sub-CC610101", "sub-CC610146", "sub-CC610292", "sub-CC620005", "sub-CC620284", - "sub-CC620413", "sub-CC620490", "sub-CC620515", "sub-CC621199", "sub-CC710037", - "sub-CC710214", "sub-CC720103", "sub-CC721392", "sub-CC721648", "sub-CC721888", - "sub-CC722421", "sub-CC722536", "sub-CC720329", - ]: - # Remove headshape points 1cm above lpa - keep = hs[2] > (lpa[2] + 10) - hs = hs[:, keep] - elif subject in [ - "sub-CC321428", "sub-CC410097", "sub-CC510076", "sub-CC510220", "sub-CC520560", - "sub-CC520597", "sub-CC520673", "sub-CC610285", "sub-CC610469", "sub-CC620429", - "sub-CC620451", "sub-CC620821", "sub-CC710494", "sub-CC722651", "sub-CC110101", - "sub-CC122172", - ]: - # Remove headshape points 2cm above lpa - keep = hs[2] > (lpa[2] + 20) - hs = hs[:, keep] - elif subject in ["sub-CC412004", "sub-CC721704"]: - # Remove headshape points 3cm above rpa - keep = hs[2] > (rpa[2] + 30) - hs = hs[:, keep] - elif subject in [ - "sub-CC110033", "sub-CC510163", "sub-CC520287", "sub-CC520607", "sub-CC620567", - ]: - # Remove headshape points 2cm below lpa - keep = hs[2] > (lpa[2] - 20) - hs = hs[:, keep] - - # Overwrite headshape file - utils.logger.log_or_print(f"overwritting {filenames['polhemus_headshape_file']}") - np.savetxt(filenames["polhemus_headshape_file"], hs) - -def custom_coregister( - src_dir, - subject, - preproc_file, - smri_file, - epoch_file, -): - if subject in [ - "sub-CC110101", "sub-CC110411", "sub-CC120049", "sub-CC420322", "sub-CC120469", - "sub-CC120208", "sub-CC120309", "sub-CC120319", "sub-CC120347", "sub-CC120376", - "sub-CC120462", "sub-CC120640", "sub-CC120727", "sub-CC120764", "sub-CC420566", - "sub-CC122405", "sub-CC410121", "sub-CC410387", "sub-CC122620", "sub-CC210023", - "sub-CC210148", "sub-CC210172", "sub-CC420888", "sub-CC210519", "sub-CC210617", - "sub-CC212153", "sub-CC220107", "sub-CC220115", "sub-CC220198", "sub-CC220234", - "sub-CC220323", "sub-CC220372", "sub-CC220567", "sub-CC220697", "sub-CC220713", - "sub-CC221002", "sub-CC221031", "sub-CC221033", "sub-CC221040", "sub-CC221107", - "sub-CC221324", "sub-CC221336", "sub-CC221352", "sub-CC221373", "sub-CC221511", - "sub-CC221828", "sub-CC221977", "sub-CC222555", "sub-CC510226", "sub-CC222652", - "sub-CC223085", "sub-CC223115", "sub-CC310008", "sub-CC310051", "sub-CC310052", - "sub-CC310214", "sub-CC310414", "sub-CC310463", "sub-CC320002", "sub-CC320059", - "sub-CC510243", "sub-CC320109", "sub-CC320202", "sub-CC320206", "sub-CC320218", - "sub-CC320325", "sub-CC320379", "sub-CC320417", "sub-CC320429", "sub-CC320500", - "sub-CC320553", "sub-CC320575", "sub-CC320576", "sub-CC320621", "sub-CC320661", - "sub-CC320686", "sub-CC320814", "sub-CC321025", "sub-CC321053", "sub-CC321107", - "sub-CC321137", "sub-CC321281", "sub-CC321291", "sub-CC321331", "sub-CC321368", - "sub-CC321368", "sub-CC321368", "sub-CC321504", "sub-CC321506", "sub-CC321529", - "sub-CC321544", "sub-CC321595", "sub-CC321880", "sub-CC321899", "sub-CC410040", - "sub-CC410091", "sub-CC410113", "sub-CC410129", "sub-CC410177", "sub-CC410289", - "sub-CC410297", "sub-CC410390", "sub-CC412004", "sub-CC420071", "sub-CC420075", - "sub-CC420143", "sub-CC420148", "sub-CC420173", "sub-CC420182", "sub-CC420226", - "sub-CC420324", "sub-CC420623", "sub-CC420776", "sub-CC321431", "sub-CC510392", - "sub-CC510393", "sub-CC510438", "sub-CC510534", "sub-CC520011", "sub-CC520134", - "sub-CC520239", "sub-CC520503", "sub-CC520517", "sub-CC520552", "sub-CC520560", - "sub-CC520673", "sub-CC520868", "sub-CC610028", "sub-CC610076", "sub-CC610594", - "sub-CC610671", "sub-CC620026", "sub-CC620073", "sub-CC620359", "sub-CC620436", - "sub-CC620935", "sub-CC621184", "sub-CC621284", "sub-CC710154", "sub-CC711035", - "sub-CC711158", "sub-CC720400", "sub-CC721291", "sub-CC721377", "sub-CC721434", - "sub-CC721707", "sub-CC121106", "sub-CC220284", "sub-CC220843", "sub-CC221740", - "sub-CC720407", - ]: - include_nose = True - use_nose = True - use_headshape = True - elif subject in [ - "sub-CC120264", "sub-CC220098", "sub-CC220999", "sub-CC410169", "sub-CC110056", - "sub-CC110174", "sub-CC112141", "sub-CC120065", "sub-CC210051", "sub-CC221565", - "sub-CC310203", - ]: - include_nose = False - use_nose = False - use_headshape = False - else: - include_nose = False - use_nose = False - use_headshape = True - - if subject in [ - "sub-CC110187", "sub-CC120276", "sub-CC310473", "sub-CC320759", "sub-CC321976", - "sub-CC410084", "sub-CC410121", "sub-CC420004", "sub-CC420094", "sub-CC420462", - "sub-CC420589", "sub-CC510050", "sub-CC510076", "sub-CC510208", "sub-CC510321", - "sub-CC510480", "sub-CC520215", "sub-CC520597", "sub-CC520624", "sub-CC610178", - "sub-CC610344", "sub-CC610392", "sub-CC610508", "sub-CC620152", "sub-CC620259", - "sub-CC620284", "sub-CC620685", "sub-CC710088", "sub-CC710429", "sub-CC710664", - "sub-CC710679", "sub-CC711244", "sub-CC720119", "sub-CC720329", - ]: - # These subjects only use a few initialisations which leads to high - # run-to-run variability. This script may need to be re-run a few times - # for these subjects to get a good coregistration. - n_init = 10 - else: - n_init = 500 - - # Compute surfaces - source_recon.wrappers.compute_surfaces( - src_dir, subject, preproc_file, smri_file, epoch_file, include_nose - ) - - # Coregister - source_recon.wrappers.coregister( - src_dir, - subject, - preproc_file, - smri_file, - epoch_file, - use_nose, - use_headshape, - n_init=n_init, - ) - -if __name__ == "__main__": - utils.logger.set_up(level="INFO") - source_recon.setup_fsl(FSL_DIR) - - # Get subjects - subjects = [] - for subject in sorted( - glob( - PREPROC_DIR - + "/mf2pt2_*_ses-rest_task-rest_meg" - + "/mf2pt2_sub-*_ses-rest_task-rest_meg_preproc_raw.fif" - ) - ): - subjects.append(pathlib.Path(subject).stem.split("_")[1]) - - # Setup files - smri_files = [] - preproc_files = [] - for subject in subjects: - smri_files.append(SMRI_FILE.format(subject)) - preproc_files.append(PREPROC_FILE.format(subject)) - - # Setup parallel processing - client = Client(n_workers=16, threads_per_worker=1) - - # Run coregistration - source_recon.run_src_batch( - config, - src_dir=COREG_DIR, - subjects=subjects, - preproc_files=preproc_files, - smri_files=smri_files, - extra_funcs=[fix_headshape_points, custom_coregister], - dask_client=True, - ) diff --git a/examples/camcan/evaluate_preproc.py b/examples/camcan/evaluate_preproc.py index 0a6a6f62..3330511a 100644 --- a/examples/camcan/evaluate_preproc.py +++ b/examples/camcan/evaluate_preproc.py @@ -2,6 +2,8 @@ """ +# Authors: Chetan Gohil + import numpy as np from glob import glob diff --git a/examples/mrc_meguk/aston/1_preprocess.py b/examples/mrc_meguk/aston/1_preprocess.py new file mode 100644 index 00000000..3341378d --- /dev/null +++ b/examples/mrc_meguk/aston/1_preprocess.py @@ -0,0 +1,39 @@ +"""Preprocesses data from the Aston site in the MRC MEGUK Partnership dataset. + +Note: there is no EOG/ECG from this site. +""" + +# Authors : Chetan Gohil + +from osl.preprocessing import run_proc_batch +from glob import glob + +# Settings +config = """ + meta: + event_codes: + visual: 30 + motor_short: 31 + motor_long: 32 + preproc: + - crop: {tmin: 10} + - find_events: {min_duration: 0.005} + - filter: {l_freq: 0.1, h_freq: 175} + - notch_filter: {freqs: 50 100 150} + - bad_channels: {picks: 'mag'} + - bad_channels: {picks: 'grad'} + - bad_segments: {segment_len: 2000, picks: 'mag'} + - bad_segments: {segment_len: 2000, picks: 'grad'} + - resample: {sfreq: 400, npad: 'auto'} +""" + +# Raw data filenames +raw_data_files = [] +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Aston/derivatives/*/meg/*resteyesclosed*.fif")) +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Aston/derivatives/*/meg/*resteyesopen*/*.fif")) +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Aston/derivatives/*/meg/*visuomotor*/*.fif")) + +# Directory for preprocessed data +output_dir = "/ohba/pi/mwoolrich/cgohil/mrc_meguk_public/Aston" + +run_proc_batch(config, raw_data_files, output_dir, nprocesses=4) diff --git a/examples/mrc_meguk/cambridge/1_preprocess.py b/examples/mrc_meguk/cambridge/1_preprocess.py new file mode 100644 index 00000000..adcee0b7 --- /dev/null +++ b/examples/mrc_meguk/cambridge/1_preprocess.py @@ -0,0 +1,49 @@ +"""Preprocessing. + +""" +from glob import glob +from pathlib import Path +from dask.distributed import Client + +from osl import preprocessing, utils + +# Authors : Rukuang Huang +# Chetan Gohil + +TASK = "resteyesopen" # resteyesopen or resteyesclosed + +RAW_DIR = "/well/woolrich/projects/mrc_meguk/raw/Cambridge/derivatives" +RAW_FILE = RAW_DIR + "/{0}/meg/{0}_task-{1}_proc-sss_meg.fif" +PREPROC_DIR = f"/well/woolrich/projects/mrc_meguk/cambridge/{TASK}/preproc" + +config = """ + preproc: + - pick: {picks: [meg, eog, ecg]} + - filter: {l_freq: 0.5, h_freq: 125, method: iir, iir_params: {order: 5, ftype: butter}} + - notch_filter: {freqs: 50 100} + - resample: {sfreq: 250} + - bad_segments: {segment_len: 500, picks: meg, significance_level: 0.1} + - bad_segments: {segment_len: 500, picks: meg, mode: diff, significance_level: 0.1} + - bad_channels: {picks: meg, significance_level: 0.1} + - ica_raw: {picks: meg, n_components: 64} + - ica_autoreject: {picks: meg, ecgmethod: correlation, eogthreshold: auto} + - interpolate_bads: {} +""" +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + client = Client(n_workers=16, threads_per_worker=1) + + inputs = [] + for directory in sorted(glob(RAW_DIR + "/sub-*")): + subject = Path(directory).name + raw_file = RAW_FILE.format(subject, TASK) + if Path(raw_file).exists(): + inputs.append(raw_file) + + preprocessing.run_proc_batch( + config, + inputs, + outdir=PREPROC_DIR, + overwrite=True, + dask_client=True, + ) diff --git a/examples/mrc_meguk/cambridge/2_coregister.py b/examples/mrc_meguk/cambridge/2_coregister.py new file mode 100644 index 00000000..2649164e --- /dev/null +++ b/examples/mrc_meguk/cambridge/2_coregister.py @@ -0,0 +1,96 @@ +"""Coregisteration. + +""" +from glob import glob +from pathlib import Path +import numpy as np +from dask.distributed import Client + +from osl import source_recon, utils + +# Authors : Rukuang Huang +# Chetan Gohil + +TASK = "resteyesopen" # resteyesopen or resteyesclosed + +BASE_DIR = "/well/woolrich/projects/mrc_meguk/cambridge/ec" +PREPROC_DIR = BASE_DIR + "/preproc" +SRC_DIR = BASE_DIR + "/src" +PREPROC_FILE = PREPROC_DIR + "/{0}_task-{1}_proc-sss_meg/{0}_task-{1}_proc-sss_meg_preproc_raw.fif" +SMRI_FILE = "/well/woolrich/projects/mrc_meguk/raw/Cambridge/{0}/anat/{0}_T1w.nii.gz" +FSL_DIR = "/well/woolrich/projects/software/fsl" + + +def fix_headshape_points(src_dir, subject, preproc_file, smri_file, epoch_file): + filenames = source_recon.rhino.get_coreg_filenames(src_dir, subject) + + # Load saved headshape and nasion files + hs = np.loadtxt(filenames["polhemus_headshape_file"]) + nas = np.loadtxt(filenames["polhemus_nasion_file"]) + lpa = np.loadtxt(filenames["polhemus_lpa_file"]) + rpa = np.loadtxt(filenames["polhemus_rpa_file"]) + + # Drop nasion by 4cm + nas[2] -= 40 + distances = np.sqrt( + (nas[0] - hs[0]) ** 2 + (nas[1] - hs[1]) ** 2 + (nas[2] - hs[2]) ** 2 + ) + keep = distances > 70 + hs = hs[:, keep] + + # Remove anything outside of rpa + keep = hs[0] < rpa[0] + hs = hs[:, keep] + + # Remove anything outside of lpa + keep = hs[0] > lpa[0] + hs = hs[:, keep] + + # if subject == "sub-cam056": + # # Remove headshape points that are 1 cm above rpa + # keep = hs[2] > (rpa[2] + 10) + # hs = hs[:, keep] + + # Overwrite headshape file + utils.logger.log_or_print(f"overwritting {filenames['polhemus_headshape_file']}") + np.savetxt(filenames["polhemus_headshape_file"], hs) + + +config = """ + source_recon: + - extract_fiducials_from_fif: {} + - fix_headshape_points: {} + - compute_surfaces_coregister_and_forward_model: + include_nose: false + use_nose: false + use_headshape: true + model: Single Layer +""" + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + source_recon.setup_fsl(FSL_DIR) + + subjects = [] + smri_files = [] + preproc_files = [] + + for directory in sorted( + glob(PREPROC_DIR + f"/sub*_task-{TASK}_proc-sss_meg") + ): + subject = Path(directory).name.split("_")[0] + subjects.append(subject) + smri_files.append(SMRI_FILE.format(subject)) + preproc_files.append(PREPROC_FILE.format(subject, TASK)) + + client = Client(n_workers=16, threads_per_worker=1) + + source_recon.run_src_batch( + config, + src_dir=SRC_DIR, + subjects=subjects, + preproc_files=preproc_files, + smri_files=smri_files, + extra_funcs=[fix_headshape_points], + dask_client=True, + ) diff --git a/examples/mrc_meguk/cambridge/3_source_reconstruct.py b/examples/mrc_meguk/cambridge/3_source_reconstruct.py new file mode 100644 index 00000000..23924ccf --- /dev/null +++ b/examples/mrc_meguk/cambridge/3_source_reconstruct.py @@ -0,0 +1,59 @@ +"""Source reconstruction (beamforming and parcellation). + +""" + +from pathlib import Path +from glob import glob +from dask.distributed import Client + +from osl import source_recon, utils + +# Authors : Rukuang Huang +# Chetan Gohil + +TASK = "resteyesopen" # resteyesopen or resteyesclosed + +BASE_DIR = "/well/woolrich/projects/mrc_meguk/cambridge/ec" +PREPROC_DIR = BASE_DIR + "/preproc" +SRC_DIR = BASE_DIR + "/src" +PREPROC_FILE = PREPROC_DIR + "/{0}_task-{1}_proc-sss_meg/{0}_task-{1}_proc-sss_meg_preproc_raw.fif" +SMRI_FILE = "/well/woolrich/projects/mrc_meguk/raw/Cambridge/{0}/anat/{0}_T1w.nii.gz" +FSL_DIR = "/well/woolrich/projects/software/fsl" + +config = """ + source_recon: + - beamform_and_parcellate: + freq_range: [1, 45] + chantypes: [mag, grad] + rank: {meg: 60} + parcellation_file: fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz + method: spatial_basis + orthogonalisation: symmetric +""" + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + source_recon.setup_fsl(FSL_DIR) + + subjects = [] + smri_files = [] + preproc_files = [] + + for directory in sorted( + glob(PREPROC_DIR + f"/sub*_task-{TASK}_proc-sss_meg") + ): + subject = Path(directory).name.split("_")[0] + subjects.append(subject) + smri_files.append(SMRI_FILE.format(subject)) + preproc_files.append(PREPROC_FILE.format(subject, task)) + + client = Client(n_workers=16, threads_per_worker=1) + + source_recon.run_src_batch( + config, + src_dir=SRC_DIR, + subjects=subjects, + preproc_files=preproc_files, + smri_files=smri_files, + dask_client=True, + ) diff --git a/examples/mrc_meguk/cambridge/4_sign_flip.py b/examples/mrc_meguk/cambridge/4_sign_flip.py new file mode 100644 index 00000000..8ee6c6de --- /dev/null +++ b/examples/mrc_meguk/cambridge/4_sign_flip.py @@ -0,0 +1,54 @@ +"""Dipole sign flipping. + +""" + +from glob import glob +from pathlib import Path +from dask.distributed import Client +from osl import utils + +from osl.source_recon import find_template_subject, run_src_batch, setup_fsl + +# Authors : Rukuang Huang +# Chetan Gohil + +TASK = "resteyesopen" # resteyesopen or resteyesclosed + +BASE_DIR = f"/well/woolrich/projects/mrc_meguk/cambridge/{TASK}" +PREPROC_DIR = BASE_DIR + "/preproc" +SRC_DIR = BASE_DIR + "/src" +FSL_DIR = "/well/woolrich/projects/software/fsl" + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + setup_fsl(FSL_DIR) + + subjects = [] + for directory in sorted( + glob(PREPROC_DIR + f"/sub*_task-{TASK}_proc-sss_meg") + ): + subject = Path(directory).name.split("_")[0] + subjects.append(subject) + + # Find a good template subject to align other subjects to + template = find_template_subject( + SRC_DIR, subjects, n_embeddings=15, standardize=True + ) + + # Settings for batch processing + config = f""" + source_recon: + - fix_sign_ambiguity: + template: {template} + n_embeddings: 15 + standardize: True + n_init: 3 + n_iter: 3000 + max_flips: 20 + """ + + # Setup parallel processing + client = Client(n_workers=16, threads_per_worker=1) + + # Do the sign flipping + run_src_batch(config, SRC_DIR, subjects, dask_client=True) diff --git a/examples/mrc_meguk/cardiff/1_preprocess.py b/examples/mrc_meguk/cardiff/1_preprocess.py new file mode 100644 index 00000000..55df0ad2 --- /dev/null +++ b/examples/mrc_meguk/cardiff/1_preprocess.py @@ -0,0 +1,40 @@ +"""Preprocesses data from the Cardiff site in the MRC MEGUK Partnership dataset. + +""" + +from osl.preprocessing import run_proc_batch +from glob import glob + +# Authors : Chetan Gohil + +# Settings +config = """ + meta: + event_codes: + visual: 30 + motor_short: 31 + motor_long: 32 + preproc: + - crop: {tmin: 10} + - set_channel_types: {EEG057: eog, EEG058: eog, EEG059: ecg} + - pick_types: {meg: true, eeg: false, eog: true, ecg: true, stim: true, ref_meg: false} + - find_events: {min_duration: 0.005} + - filter: {l_freq: 0.1, h_freq: 175} + - notch_filter: {freqs: 50 100 150} + - bad_channels: {picks: 'meg'} + - bad_segments: {segment_len: 2000, picks: 'meg'} + - resample: {sfreq: 400, npad: 'auto'} + - ica_raw: {n_components: 60, picks: data} + - ica_autoreject: {picks: 'mag', ecgmethod: 'correlation', apply: false} +""" + +# Raw data filenames +raw_data_files = [] +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Cardiff/*/meg/*resteyesclosed*/*.meg4")) +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Cardiff/*/meg/*resteyesopen*/*.meg4")) +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Cardiff/*/meg/*visuomotor*/*.meg4")) + +# Directory for preprocessed data +output_dir = "/ohba/pi/mwoolrich/cgohil/mrc_meguk_public/Cardiff" + +run_proc_batch(config, raw_data_files, output_dir, nprocesses=4) diff --git a/examples/mrc_meguk/glasgow/1_preprocess.py b/examples/mrc_meguk/glasgow/1_preprocess.py new file mode 100644 index 00000000..2005edd7 --- /dev/null +++ b/examples/mrc_meguk/glasgow/1_preprocess.py @@ -0,0 +1,37 @@ +"""Preprocesses data from the Glasgow site in the MRC MEGUK Partnership dataset. + +""" + +from osl.preprocessing import run_proc_batch +from glob import glob + +# Authors : Chetan Gohil + +# Settings +config = """ + meta: + event_codes: + visual: 30 + motor_short: 31 + motor_long: 32 + preproc: + - crop: {tmin: 10} + - find_events: {min_duration: 0.005} + - filter: {l_freq: 0.1, h_freq: 175} + - notch_filter: {freqs: 50 100 150} + - bad_channels: {picks: 'mag'} + - bad_segments: {segment_len: 2000, picks: 'mag'} + - resample: {sfreq: 400, npad: 'auto'} + - ica_raw: {n_components: 60, picks: 'meg'} +""" + +# Raw data filenames +raw_data_files = [] +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Glasgow/*/meg/*resteyesclosed*/c,rfDC")) +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Glasgow/*/meg/*resteyesopen*/c,rfDC")) +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Glasgow/*/meg/*visuomotor*/c,rfDC")) + +# Directory for preprocessed data +output_dir = "/ohba/pi/mwoolrich/cgohil/mrc_meguk_public/Glasgow" + +run_proc_batch(config, raw_data_files, output_dir, nprocesses=4) diff --git a/examples/mrc_meguk/notts/1_preprocess.py b/examples/mrc_meguk/notts/1_preprocess.py new file mode 100755 index 00000000..daac66b3 --- /dev/null +++ b/examples/mrc_meguk/notts/1_preprocess.py @@ -0,0 +1,52 @@ +"""Preprocessing. + +""" + +from glob import glob +from pathlib import Path +from dask.distributed import Client + +from osl import preprocessing, utils + +# Authors : Rukuang Huang +# Chetan Gohil + +TASK = "resteyesopen" # resteyesopen or resteyesclosed + +RAW_DIR = "/well/woolrich/projects/mrc_meguk/raw/Nottingham" +RAW_FILE = RAW_DIR + "/{0}/meg/{0}_task-{1}_meg.ds" +PREPROC_DIR = f"/well/woolrich/projects/mrc_meguk/notts/{TASK}/preproc" + +config = """ + preproc: + - set_channel_types: {EEG057: eog, EEG058: eog, EEG059: ecg} + - pick: {picks: [mag, eog, ecg]} + - filter: {l_freq: 0.5, h_freq: 125, method: iir, iir_params: {order: 5, ftype: butter}} + - notch_filter: {freqs: 50 100} + - resample: {sfreq: 250} + - bad_segments: {segment_len: 500, picks: mag, significance_level: 0.1} + - bad_segments: {segment_len: 500, picks: mag, mode: diff, significance_level: 0.1} + - bad_channels: {picks: mag, significance_level: 0.1} + - ica_raw: {picks: mag, n_components: 64} + - ica_autoreject: {picks: mag, ecgmethod: correlation, eogthreshold: auto} + - interpolate_bads: {} +""" + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + client = Client(n_workers=16, threads_per_worker=1) + + inputs = [] + for directory in sorted(glob(RAW_DIR + "/sub-*")): + subject = Path(directory).name + raw_file = RAW_FILE.format(subject, TASK) + if Path(raw_file).exists(): + inputs.append(raw_file) + + preprocessing.run_proc_batch( + config, + inputs, + outdir=PREPROC_DIR, + overwrite=True, + dask_client=True, + ) diff --git a/examples/mrc_meguk/notts/fix_smri_files.py b/examples/mrc_meguk/notts/2_fix_smri_files.py old mode 100644 new mode 100755 similarity index 67% rename from examples/mrc_meguk/notts/fix_smri_files.py rename to examples/mrc_meguk/notts/2_fix_smri_files.py index b34b7afa..7110779c --- a/examples/mrc_meguk/notts/fix_smri_files.py +++ b/examples/mrc_meguk/notts/2_fix_smri_files.py @@ -1,5 +1,4 @@ -"""The sform headers in the structural MRIs in the MRC MEG UK dataset have an issue. -This script writes new MRI files compatible with oslpy. +"""Fix structural MRI files. """ @@ -13,22 +12,25 @@ from osl import source_recon +# Authors : Rukuang Huang +# Chetan Gohil -PREPROC_DIR = "/ohba/pi/mwoolrich/cgohil/ukmp_notts/preproc" -SMRI_FILE = ( - "/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Nottingham/{0}/anat/{0}_T1w.nii.gz" -) +TASK = "resteyesopen" # resteyesopen or resteyesclosed -FIXED_SMRI_DIR = "/ohba/pi/mwoolrich/cgohil/ukmp_notts/smri" +PREPROC_DIR = f"/well/woolrich/projects/mrc_meguk/notts/{TASK}/preproc" +SMRI_FILE = "/well/woolrich/projects/mrc_meguk/raw/Nottingham/{0}/anat/{0}_T1w.nii.gz" +FIXED_SMRI_DIR = f"/well/woolrich/projects/mrc_meguk/notts/{TASK}/smri" # Setup FSL -source_recon.setup_fsl("/home/cgohil/local/fsl") +source_recon.setup_fsl("/well/woolrich/projects/software/fsl") # Look up which subjects we preprocessed to see what SMRI files we need to fix smri_files = [] for path in sorted(glob(PREPROC_DIR + "/*/sub-*_preproc_raw.fif")): subject = Path(path).stem.split("_")[0] - smri_files.append(SMRI_FILE.format(subject)) + smri_file = SMRI_FILE.format(subject) + if Path(smri_file).exists(): + smri_files.append(smri_file) # Make output directory if it doesn't exist makedirs(FIXED_SMRI_DIR, exist_ok=True) diff --git a/examples/mrc_meguk/notts/3_source_reconstruct.py b/examples/mrc_meguk/notts/3_source_reconstruct.py new file mode 100755 index 00000000..dd67bfee --- /dev/null +++ b/examples/mrc_meguk/notts/3_source_reconstruct.py @@ -0,0 +1,126 @@ +"""Source reconstruction. + +""" + +import os.path as op +from glob import glob +from pathlib import Path + +import numpy as np +import pandas as pd +from dask.distributed import Client + +from osl import source_recon, utils + +# Authors : Rukuang Huang +# Chetan Gohil + +TASK = "resteyesopen" # resteyesopen or resteyesclosed + +# Directories +RAW_DIR = "/well/woolrich/projects/mrc_meguk/raw/Nottingham" +PREPROC_DIR = f"/well/woolrich/projects/mrc_meguk/notts/{TASK}/preproc" +SRC_DIR = f"/well/woolrich/projects/mrc_meguk/notts/{TASK}/src" + +SMRI_FILE = f"/well/woolrich/projects/mrc_meguk/notts/{TASK}/smri/{0}_T1w.nii.gz" +PREPROC_FILE = PREPROC_DIR + "/{0}_task-resteyesopen_meg/{0}_task-" + TASK + "_meg_preproc_raw.fif" +POS_FILE = RAW_DIR + "/{0}/meg/{0}_headshape.pos" + +# Setup FSL +source_recon.setup_fsl("/well/woolrich/projects/software/fsl") + +def save_polhemus_from_pos(src_dir, subject, preproc_file, smri_file, epoch_file): + """Saves fiducials/headshape from a pos file.""" + + # Load pos file + pos_file = POS_FILE.format(subject) + utils.logger.log_or_print(f"Saving polhemus from {pos_file}") + + # Get coreg filenames + filenames = source_recon.rhino.get_coreg_filenames(src_dir, subject) + + # Load in txt file, these values are in cm in polhemus space: + num_headshape_pnts = int(pd.read_csv(pos_file, header=None).to_numpy()[0]) + data = pd.read_csv(pos_file, header=None, skiprows=[0], delim_whitespace=True) + + # RHINO is going to work with distances in mm + # So convert to mm from cm, note that these are in polhemus space + data.iloc[:, 1:4] = data.iloc[:, 1:4] * 10 + + # Polhemus fiducial points in polhemus space + polhemus_nasion = ( + data[data.iloc[:, 0].str.match("nasion")] + .iloc[0, 1:4] + .to_numpy() + .astype("float64") + .T + ) + polhemus_rpa = ( + data[data.iloc[:, 0].str.match("right")] + .iloc[0, 1:4] + .to_numpy() + .astype("float64") + .T + ) + polhemus_lpa = ( + data[data.iloc[:, 0].str.match("left")] + .iloc[0, 1:4] + .to_numpy() + .astype("float64") + .T + ) + + # Polhemus headshape points in polhemus space in mm + polhemus_headshape = ( + data[0:num_headshape_pnts].iloc[:, 1:4].to_numpy().astype("float64").T + ) + + # Save + np.savetxt(filenames["polhemus_nasion_file"], polhemus_nasion) + np.savetxt(filenames["polhemus_rpa_file"], polhemus_rpa) + np.savetxt(filenames["polhemus_lpa_file"], polhemus_lpa) + np.savetxt(filenames["polhemus_headshape_file"], polhemus_headshape) + + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + client = Client(n_workers=16, threads_per_worker=1) + + # Settings + config = """ + source_recon: + - save_polhemus_from_pos: {} + - compute_surfaces_coregister_and_forward_model: + include_nose: true + use_nose: true + use_headshape: true + model: Single Layer + - beamform_and_parcellate: + freq_range: [1, 45] + chantypes: mag + rank: {mag: 120} + parcellation_file: fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz + method: spatial_basis + orthogonalisation: symmetric + """ + + # Get input files + subjects = [] + smri_files = [] + preproc_files = [] + for path in sorted(glob(PREPROC_FILE.replace("{0}", "*"))): + subject = Path(path).stem.split("_")[0] + subjects.append(subject) + preproc_files.append(PREPROC_FILE.format(subject)) + smri_files.append(SMRI_FILE.format(subject)) + + # Source reconstruction + source_recon.run_src_batch( + config, + src_dir=SRC_DIR, + subjects=subjects, + preproc_files=preproc_files, + smri_files=smri_files, + extra_funcs=[save_polhemus_from_pos], + dask_client=True, + ) diff --git a/examples/mrc_meguk/notts/4_sign_flip.py b/examples/mrc_meguk/notts/4_sign_flip.py new file mode 100755 index 00000000..1a3700eb --- /dev/null +++ b/examples/mrc_meguk/notts/4_sign_flip.py @@ -0,0 +1,49 @@ +"""Sign flipping. + +""" + +from glob import glob +from dask.distributed import Client + +from osl import source_recon, utils + +# Authors : Rukuang Huang +# Chetan Gohil + +TASK = "resteyesopen" # resteyesopen or resteyesclosed + +# Setup FSL +source_recon.setup_fsl("/well/woolrich/projects/software/fsl") + +# Directories +SRC_DIR = f"/well/woolrich/projects/mrc_meguk/notts/{TASK}/src" + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + client = Client(n_workers=16, threads_per_worker=1) + + # Subjects to sign flip + subjects = [] + for path in sorted(glob(SRC_DIR + "/*/parc/parc-raw.fif")): + subject = path.split("/")[-3] + subjects.append(subject) + + # Find a good template subject to align other subjects to + template = source_recon.find_template_subject( + SRC_DIR, subjects, n_embeddings=15, standardize=True + ) + + # Settings + config = f""" + source_recon: + - fix_sign_ambiguity: + template: {template} + n_embeddings: 15 + standardize: True + n_init: 3 + n_iter: 3000 + max_flips: 20 + """ + + # Do the sign flipping + source_recon.run_src_batch(config, SRC_DIR, subjects, dask_client=True) diff --git a/examples/mrc_meguk/notts/preproc_and_parcellate.py b/examples/mrc_meguk/notts/preproc_and_parcellate.py deleted file mode 100644 index c2c9298f..00000000 --- a/examples/mrc_meguk/notts/preproc_and_parcellate.py +++ /dev/null @@ -1,254 +0,0 @@ -#!/usr/bin/env python - -"""Full pipeline for getting parcel time series from Nottingham MEG UK partnership data. - -Includes preprocessing, beamforming and parcellation. -""" - -# Authors: Mark Woolrich -# Chetan Gohil - -import os -import os.path as op - -import numpy as np -import pandas as pd - -from osl import preprocessing -from osl import source_recon - -subjects_dir = "/Users/woolrich/homedir/vols_data/ukmp" -subjects_to_do = np.arange(75)+1 -subjects_to_do = np.arange(2)+1 - -task = "resteyesopen" - -subjects = [] -ds_files = [] -preproc_fif_files = [] -smri_files = [] - -recon_dir = op.join(subjects_dir, 'recon') -POS_FILE = subjects_dir + "/{0}/meg/{0}_headshape.pos" - -# input files -for sub in subjects_to_do: - subject = "sub-not" + ("{}".format(sub)).zfill(3) - ds_file = op.join( - subjects_dir, subject, "meg", subject + "_task-" + task + "_meg.ds" - ) - smri_file = op.join(subjects_dir, subject, "anat", subject + "_T1w.nii.gz") - - # check ds file and structural file exists for this subject - if op.exists(ds_file) and op.exists(smri_file): - ds_files.append(ds_file) - subjects.append(subject) - preproc_fif_files.append( - op.join( - subjects_dir, - subject, - "meg", - subject + "_task-" + task + "_meg_preproc_raw.fif", - ) - ) - - smri_files.append(smri_file) - -run_preproc = True -run_beamform_and_parcellate = True -run_fix_sign_ambiguity = False - -# parcellation to use -parcellation_fname = "fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz" - -# resolution of dipole grid for source recon -gridstep = 7 # mm - -if run_preproc: - - # Note that with CTF mne.pick_types will return: - # ~274 axial grads (as magnetometers) if {picks: 'mag', ref_meg: False} - # ~28 reference axial grads if {picks: 'grad'} - - config = """ - preproc: - - resample: {sfreq: 250} - - filter: {l_freq: 0.5, h_freq: 100, method: 'iir', iir_params: {order: 5, ftype: butter}} - - crop: {tmin: 10, tmax: 300} - - set_channel_types: {EEG057: eog, EEG058: eog, EEG059: ecg} - - bad_channels: {picks: 'mag', ref_meg: False, significance_level: 0.1} - - bad_channels: {picks: 'grad', significance_level: 0.4} - - bad_segments: {segment_len: 600, picks: 'mag', ref_meg: False, significance_level: 0.1} - """ - - # - crop: {tmin: 20, tmax: 280} - #- bad_channels: {picks: 'meg', significance_level: 0.4} - #- bad_segments: {segment_len: 600, picks: 'meg', significance_level: 0.1} - - # Process a single file, this outputs fif_file - dataset = preprocessing.run_proc_batch( - config, ds_files, outdir=subjects_dir, overwrite=True - ) - - for subject in subjects: - os.system( - "mv {} {}".format( - op.join(subjects_dir, subject + "_task-" + task + "_meg_preproc_raw.*"), - op.join(subjects_dir, subject, "meg"), - ) - ) - -# ------------------------------------------------------------- -# %% Coreg and Source recon and Parcellate - -def save_polhemus_from_pos(src_dir, subject, preproc_file, smri_file, tmp): - """Saves fiducials/headshape from a pos file.""" - - # Load pos file - pos_file = POS_FILE.format(subject) - - #  Get coreg filenames - filenames = source_recon.rhino.get_coreg_filenames(src_dir, subject) - - # Load in txt file, these values are in cm in polhemus space: - num_headshape_pnts = int(pd.read_csv(pos_file, header=None).to_numpy()[0]) - data = pd.read_csv(pos_file, header=None, skiprows=[0], delim_whitespace=True) - - # RHINO is going to work with distances in mm - # So convert to mm from cm, note that these are in polhemus space - data.iloc[:, 1:4] = data.iloc[:, 1:4] * 10 - - # Polhemus fiducial points in polhemus space - polhemus_nasion = ( - data[data.iloc[:, 0].str.match("nasion")] - .iloc[0, 1:4].to_numpy().astype("float64").T - ) - polhemus_rpa = ( - data[data.iloc[:, 0].str.match("right")] - .iloc[0, 1:4].to_numpy().astype("float64").T - ) - polhemus_lpa = ( - data[data.iloc[:, 0].str.match("left")] - .iloc[0, 1:4].to_numpy().astype("float64").T - ) - - # Polhemus headshape points in polhemus space in mm - polhemus_headshape = ( - data[0:num_headshape_pnts] - .iloc[:, 1:4].to_numpy().astype("float64").T - ) - - # Save - np.savetxt(filenames["polhemus_nasion_file"], polhemus_nasion) - np.savetxt(filenames["polhemus_rpa_file"], polhemus_rpa) - np.savetxt(filenames["polhemus_lpa_file"], polhemus_lpa) - np.savetxt(filenames["polhemus_headshape_file"], polhemus_headshape) - -if run_beamform_and_parcellate: - - # Settings - config = """ - source_recon: - - save_polhemus_from_pos: {} - - compute_surfaces: - include_nose: true - - coregister: - use_nose: true - use_headshape: true - - forward_model: - model: Single Layer - - beamform_and_parcellate: - freq_range: [1, 45] - chantypes: mag - rank: {mag: 120} - parcellation_file: fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz - method: spatial_basis - orthogonalisation: symmetric - """ - - source_recon.run_src_batch( - config, - src_dir=recon_dir, - subjects=subjects, - preproc_files=preproc_fif_files, - smri_files=smri_files, - extra_funcs=[save_polhemus_from_pos], - ) - - -if False: - - # to just run surfaces for subject 0: - source_recon.rhino.compute_surfaces( - smri_files[0], - recon_dir, - subjects[0], - overwrite=True - ) - - # to view surfaces for subject 0: - source_recon.rhino.surfaces_display(recon_dir, subjects[3]) - - # to just run coreg for subject 0: - source_recon.rhino.coreg( - preproc_fif_files[0], - recon_dir, - subjects[0], - already_coregistered=True - ) - - # to view coreg result for subject 0: - source_recon.rhino.coreg_display(recon_dir, subjects[0], - plot_type='surf') - -if False: - # ------------------------------------------------------------- - # %% Take a look at leadfields - - # load forward solution - fwd_fname = rhino.get_coreg_filenames(subjects_dir, subject)['forward_model_file'] - fwd = mne.read_forward_solution(fwd_fname) - - leadfield = fwd['sol']['data'] - print("Leadfield size : %d sensors x %d dipoles" % leadfield.shape) - -# ------------------------------------------------------------- -# %% Sign flip - -if run_fix_sign_ambiguity: - # Find a good template subject to align other subjects to - template = source_recon.find_template_subject( - recon_dir, subjects, n_embeddings=15, standardize=True - ) - - # Settings for batch processing - config = f""" - source_recon: - - fix_sign_ambiguity: - template: {template} - n_embeddings: 13 - standardize: True - n_init: 3 - n_iter: 2500 - max_flips: 20 - """ - - # Do the sign flipping - source_recon.run_src_batch(config, recon_dir, subjects, report_name='sflip_report') - - if True: - # copy sf files to a single directory (makes it easier to copy minimal files to, e.g. BMRC, for downstream analysis) - os.makedirs(op.join(recon_dir, 'sflip_data'), exist_ok=True) - - sflip_parc_files = [] - for subject in subjects: - sflip_parc_file_from = op.join(recon_dir, subject, 'sflip_parc-raw.fif') - sflip_parc_file_to = op.join(recon_dir, 'sflip_data', subject + '_sflip_parc-raw.fif') - - os.system('cp -f {} {}'.format(sflip_parc_file_from, sflip_parc_file_to)) - - sflip_parc_files.append(sflip_parc_file_to) - - - - diff --git a/examples/mrc_meguk/oxford/1_preprocess.py b/examples/mrc_meguk/oxford/1_preprocess.py new file mode 100644 index 00000000..95df1c7d --- /dev/null +++ b/examples/mrc_meguk/oxford/1_preprocess.py @@ -0,0 +1,40 @@ +"""Preprocesses data from the Oxford site in the MRC MEGUK Partnership dataset. + +""" + +from osl.preprocessing import run_proc_batch +from glob import glob + +# Authors : Chetan Gohil + +# Settings +config = """ + meta: + event_codes: + visual: 30 + motor_short: 31 + motor_long: 32 + preproc: + - crop: {tmin: 10} + - find_events: {min_duration: 0.005} + - filter: {l_freq: 0.1, h_freq: 175} + - notch_filter: {freqs: 50 100 150} + - bad_channels: {picks: 'mag'} + - bad_channels: {picks: 'grad'} + - bad_segments: {segment_len: 2000, picks: 'mag'} + - bad_segments: {segment_len: 2000, picks: 'grad'} + - resample: {sfreq: 400, npad: 'auto'} + - ica_raw: {n_components: 60, picks: 'meg'} + - ica_autoreject: {picks: 'meg', ecgmethod: 'correlation', apply: false} +""" + +# Raw data filenames +raw_data_files = [] +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Oxford/*/meg/*resteyesclosed*.fif")) +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Oxford/*/meg/*resteyesopen*/*.fif")) +raw_data_files += sorted(glob("/ohba/pi/mwoolrich/datasets/mrc_meguk_public/Oxford/*/meg/*visuomotor*/*.fif")) + +# Directory for preprocessed data +output_dir = "/ohba/pi/mwoolrich/cgohil/mrc_meguk_public/Oxford" + +run_proc_batch(config, raw_data_files, output_dir, nprocesses=4) diff --git a/examples/oxford/1_maxfilter.py b/examples/oxford/1_maxfilter.py new file mode 100644 index 00000000..9f6390ec --- /dev/null +++ b/examples/oxford/1_maxfilter.py @@ -0,0 +1,28 @@ +"""Example script for maxfiltering raw data recorded at Oxford using the new scanner. + +Note: this script needs to be run on a computer with a MaxFilter license. +""" + +# Authors: Chetan Gohil + +from osl.maxfilter import run_maxfilter_batch + +# Setup paths to raw (pre-maxfiltered) fif files +input_files = [ + "path/to/file1.fif" + "path/to/file2.fif" +] + +# Directory to save the maxfiltered data to +output_directory = "output/maxfilter" + +# Run MaxFiltering +# +# Note: +# - We don't use the -trans option because it affects the placement of the head during coregistration. +# - See the /maxfilter directory for further info. +run_maxfilter_batch( + input_files, + output_directory, + "--maxpath /neuro/bin/util/maxfilter --scanner Neo --tsss --mode multistage --headpos --movecomp", +) diff --git a/examples/oxford/2_preprocess.py b/examples/oxford/2_preprocess.py new file mode 100644 index 00000000..0bb95d94 --- /dev/null +++ b/examples/oxford/2_preprocess.py @@ -0,0 +1,57 @@ +"""Preprocessing. + +Note, we preprocess multiple subjects in parallel to speed things up. +""" + +# Authors: Chetan Gohil + +from dask.distributed import Client + +from osl import preprocessing, utils + +# Files and directories +raw_file = "output/maxfilter/{subject}_tsss.fif" # {subject} will be replace by the name for the subject +preproc_dir = "output/preproc" # output directory containing the preprocess files + +subjects = ["sub-001", "sub-002"] + +# Settings +config = """ + preproc: + - crop: {tmin: 30} + - filter: {l_freq: 0.5, h_freq: 125, method: iir, iir_params: {order: 5, ftype: butter}} + - notch_filter: {freqs: 50 100} + - resample: {sfreq: 250} + - bad_segments: {segment_len: 500, picks: mag, significance_level: 0.1} + - bad_segments: {segment_len: 500, picks: grad, significance_level: 0.1} + - bad_segments: {segment_len: 500, picks: mag, mode: diff, significance_level: 0.1} + - bad_segments: {segment_len: 500, picks: grad, mode: diff, significance_level: 0.1} + - bad_channels: {picks: mag, significance_level: 0.1} + - bad_channels: {picks: grad, significance_level: 0.1} + - ica_raw: {picks: meg, n_components: 0.99} + - ica_autoreject: {picks: meg, ecgmethod: correlation, eogthreshold: auto} + - interpolate_bads: {} +""" + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + + # Get input files + inputs = [] + for subject in subjects: + inputs.append(raw_file.format(subject=subject)) + + # Setup parallel processing + # + # n_workers is the number of CPUs to use, + # we recommend less than half the total number of CPUs you have + client = Client(n_workers=4, threads_per_worker=1) + + # Main preprocessing + preprocessing.run_proc_batch( + config, + inputs, + outdir=preproc_dir, + overwrite=True, + dask_client=True, + ) diff --git a/examples/oxford/3_coregister.py b/examples/oxford/3_coregister.py new file mode 100644 index 00000000..71dd47be --- /dev/null +++ b/examples/oxford/3_coregister.py @@ -0,0 +1,102 @@ +"""Coregisteration. + +The scripts was first run for all subjects (with n_init=1). Then for subjects +whose coregistration looked a bit off we re-run this script just for that +particular subject with a higher n_init. + +Note, these scripts do not include/use the nose in the coregistration. +If you want to use the nose you need to change the config to include the nose +and you may not want to call the fix_headshape_points function. +""" + +# Authors: Chetan Gohil + +import numpy as np +from dask.distributed import Client + +from osl import source_recon, utils + +# Directories +preproc_dir = "/output/preproc" +anat_dir = "path/to/smri/dir" +coreg_dir = "output/coreg" +fsl_dir = "/opt/ohba/fsl/6.0.5" # this is where FSL is installed on hbaws + +# Files ({subject} will be replaced by the name for the subject) +preproc_file = preproc_dir + "/{subject}_tsss_preproc_raw.fif" +smri_file = anat_dir + "/{subject}/anat/{subject}_T1w.nii" + +# Subjects to coregister +subjects = ["sub-001", "sub-002"] + +# Settings +config = """ + source_recon: + - extract_fiducials_from_fif: {} + - fix_headshape_points: {} + - compute_surfaces: + include_nose: False + - coregister: + use_nose: False + use_headshape: True + #n_init: 50 +""" + +def fix_headshape_points(src_dir, subject, preproc_file, smri_file, epoch_file): + filenames = source_recon.rhino.get_coreg_filenames(src_dir, subject) + + # Load saved headshape and nasion files + hs = np.loadtxt(filenames["polhemus_headshape_file"]) + nas = np.loadtxt(filenames["polhemus_nasion_file"]) + lpa = np.loadtxt(filenames["polhemus_lpa_file"]) + rpa = np.loadtxt(filenames["polhemus_rpa_file"]) + + # Remove headshape points on the nose + remove = np.logical_and(hs[1] > max(lpa[1], rpa[1]), hs[2] < nas[2]) + hs = hs[:, ~remove] + + # Remove headshape points on the neck + remove = hs[2] < min(lpa[2], rpa[2]) - 4 + hs = hs[:, ~remove] + + # Remove headshape points far from the head in any direction + remove = np.logical_or( + hs[0] < lpa[0] - 5, + np.logical_or( + hs[0] > rpa[0] + 5, + hs[1] > nas[1] + 5, + ), + ) + hs = hs[:, ~remove] + + # Overwrite headshape file + utils.logger.log_or_print(f"overwritting {filenames['polhemus_headshape_file']}") + np.savetxt(filenames["polhemus_headshape_file"], hs) + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + source_recon.setup_fsl(fsl_dir) + + # Setup files + preproc_files = [] + smri_files = [] + for subject in subjects: + preproc_files.append(preproc_file.format(subject=subject)) + smri_files.append(smri_file.format(subject=subject)) + + # Setup parallel processing + # + # n_workers is the number of CPUs to use, + # we recommend less than half the total number of CPUs you have + client = Client(n_workers=4, threads_per_worker=1) + + # Run coregistration + source_recon.run_src_batch( + config, + src_dir=coreg_dir, + subjects=subjects, + preproc_files=preproc_files, + smri_files=smri_files, + extra_funcs=[fix_headshape_points], + dask_client=True, + ) diff --git a/examples/oxford/4_source_reconstruct.py b/examples/oxford/4_source_reconstruct.py new file mode 100644 index 00000000..45724793 --- /dev/null +++ b/examples/oxford/4_source_reconstruct.py @@ -0,0 +1,62 @@ +"""Source reconstruction. + +This includes beamforming, parcellation and orthogonalisation. + +Note, before this script is run the /coreg directory created by 3_coregister.py +must be copied and renamed to /src. +""" + +# Authors: Chetan Gohil + +from dask.distributed import Client + +from osl import source_recon, utils + +# Directories +preproc_dir = "output/preproc" +src_dir = "output/src" +fsl_dir = "/opt/ohba/fsl/6.0.5" # this is where FSL is installed on hbaws + +# Files +preproc_file = preproc_dir + "/{subject}_tsss_preproc_raw.fif" # {subject} will be replaced by the subject name + +# Subjects to do +subjects = ["sub-001", "sub-002"] + +# Settings +config = """ + source_recon: + - forward_model: + model: Single Layer + - beamform_and_parcellate: + freq_range: [1, 45] + chantypes: [mag, grad] + rank: {meg: 60} + parcellation_file: fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz + method: spatial_basis + orthogonalisation: symmetric +""" + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + source_recon.setup_fsl(fsl_dir) + + # Get paths to files + preproc_files = [] + for subject in subjects: + preproc_files.append(preproc_file.format(subject=subject)) + + # Setup parallel processing + # + # n_workers is the number of CPUs to use, + # we recommend less than half the total number of CPUs you have + client = Client(n_workers=4, threads_per_worker=1) + + # Source reconstruction + source_recon.run_src_batch( + config, + src_dir=src_dir, + subjects=subjects, + preproc_files=preproc_files, + dask_client=True, + ) diff --git a/examples/oxford/5_sign_flip.py b/examples/oxford/5_sign_flip.py new file mode 100644 index 00000000..862e9b94 --- /dev/null +++ b/examples/oxford/5_sign_flip.py @@ -0,0 +1,57 @@ +"""Performs sign flipping. + +""" + +# Authors: Chetan Gohil + +from glob import glob +from dask.distributed import Client + +from osl import utils +from osl.source_recon import find_template_subject, run_src_batch, setup_fsl + +# Directories +src_dir = "output/src" +fsl_dir = "/opt/ohba/fsl/6.0.5" # this is where FSL is installed on hbaws + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + setup_fsl(fsl_dir) + + # Subjects to sign flip + # We create a list by looking for subjects that have a parc/parc-raw.fif file + subjects = [] + for path in sorted(glob(src_dir + "/*/parc/parc-raw.fif")): + subject = path.split("/")[-3] + subjects.append(subject) + + # Find a good template subject to align other subjects to + template = find_template_subject( + src_dir, subjects, n_embeddings=15, standardize=True + ) + + # Settings for batch processing + config = f""" + source_recon: + - fix_sign_ambiguity: + template: {template} + n_embeddings: 15 + standardize: True + n_init: 3 + n_iter: 500 + max_flips: 20 + """ + + # Setup parallel processing + # + # n_workers is the number of CPUs to use, + # we recommend less than half the total number of CPUs you have + client = Client(n_workers=4, threads_per_worker=1) + + # Do the sign flipping + run_src_batch( + config, + src_dir=src_dir, + subjects=subjects + dask_client=True, + ) diff --git a/examples/oxford/README.md b/examples/oxford/README.md new file mode 100644 index 00000000..ccfc72b3 --- /dev/null +++ b/examples/oxford/README.md @@ -0,0 +1,24 @@ +Preprocessing and Source Reconstruction for Data Collected at OHBA, Oxford +-------------------------------------------------------------------------- + +The scripts in this directory contain the recommended settings for preprocessing continuous data collected using the new scanner (Neo) at OHBA, Oxford. It doesn't matter if the data was recorded during a task or at rest. The subdirectories contain scripts used to preprocess/source reconstruct data for various projects at OHBA. + +To run these scripts you need: + +1. The raw `.fif` file containing the recording for each session. + +2. A structural MRI for each subject. + +The steps are: + +1. **MaxFilter** (`1_maxfilter.py`). This script MaxFilters the raw recordings. + +2. **Preprocess** (`2_preprocess.py`). This script filters, downsamples, and uses automated algorithms to detect bad segments/channels and uses ICA to remove eye blinks. Note, the automated ICA artefact removal might not always work very well, so it's often worthwhile to check the preprocessing. + +3. **Coregistration** (`3_coregister.py`). This script aligns the sMRI and MEG space (using the polhemus head space). Here, we advise you check the head has been placed in a plausible location in the MEG scanner. You can do this by looking at the coregistration panel in the report (`coreg/report/subjects_report.html` and `summary_report.html`). You may need to re-run this script with different settings (e.g. by increasing `n_init`) to fix poorly aligned subjects. + +4. **Source reconstruct** (`4_source_reconstruct.py`). This script calculates the forward model, beamforms and parcellates the data. The parcellated data files can be found in `src/{subject}/parc/parc-raw.fif`. + +5. **Dipole sign flipping** (`5_sign_flip.py`). This script tries to minimise the effect of the sign of each parcel time course being misaligned between subjects when fitting group-level models, such as a Hidden Markov Model. This step maybe skipped if you're not interested in fitting group models. + +The final data can be found in `src/{subject}/sflip_parc-raw.fif`. diff --git a/examples/oxford/pd_gripper/README.md b/examples/oxford/pd_gripper/README.md index ee35581c..24627479 100644 --- a/examples/oxford/pd_gripper/README.md +++ b/examples/oxford/pd_gripper/README.md @@ -1,7 +1,7 @@ Example scripts to preprocess and source reconstruct task data -------------------------------------------------------------- -These scripts differ from dgresch_int_ext in the order in which we epoch the data. In these example scripts we: +These scripts differ from int_ext in the order in which we epoch the data. In these example scripts we: - Preprocess the continuous sensor data: `1_preprocess.py`. - Perform manual ICA: `2_manual_ica.py`. diff --git a/examples/wakeman_henson/preprocess.py b/examples/wakeman_henson/2023_workshop/1_preprocess.py similarity index 95% rename from examples/wakeman_henson/preprocess.py rename to examples/wakeman_henson/2023_workshop/1_preprocess.py index cde0e214..56ae95fa 100644 --- a/examples/wakeman_henson/preprocess.py +++ b/examples/wakeman_henson/2023_workshop/1_preprocess.py @@ -1,7 +1,4 @@ -#!/usr/bin/env python - -""" -Run group analysis on parcellated data on the Wakeman-Henson dataset. +"""Preprocessing. """ @@ -27,8 +24,8 @@ #subj_sess_2exclude = np.ones(subj_sess_2exclude.shape).astype(bool) #subj_sess_2exclude[0:1,0:2]=False -# ------------------------------------------------------------- -# %% Setup file names +# ---------------- +# Setup file names smri_files = [] fif_files = [] @@ -111,4 +108,3 @@ osl.preprocessing.run_proc_batch( config, fif_files, outdir=subjects_dir, overwrite=True ) - diff --git a/examples/wakeman_henson/recon_and_parcellate.py b/examples/wakeman_henson/2023_workshop/2_recon_and_parcellate.py similarity index 90% rename from examples/wakeman_henson/recon_and_parcellate.py rename to examples/wakeman_henson/2023_workshop/2_recon_and_parcellate.py index 97781fcb..5dd02f6e 100644 --- a/examples/wakeman_henson/recon_and_parcellate.py +++ b/examples/wakeman_henson/2023_workshop/2_recon_and_parcellate.py @@ -1,7 +1,4 @@ -#!/usr/bin/env python - -""" -Run group analysis on parcellated data on the Wakeman-Henson dataset. +"""Run group analysis on parcellated data on the Wakeman-Henson dataset. """ @@ -27,8 +24,8 @@ #subj_sess_2exclude = np.ones(subj_sess_2exclude.shape).astype(bool) #subj_sess_2exclude[0:1,0:1]=False -# ------------------------------------------------------------- -# %% Setup file names +# ---------------- +# Setup file names smri_files = [] preproc_fif_files = [] @@ -62,8 +59,8 @@ preproc_fif_files.append(preproc_fif_file) sflip_parc_files.append(sflip_parc_file) -# ------------------------------------------------------------- -# %% Coreg and Source recon and Parcellate +# ------------------------------------- +# Coreg and Source recon and Parcellate config = """ source_recon: @@ -94,8 +91,8 @@ smri_files=smri_files, ) -# ------------------------------------------------------------- -# %% Sign flip +# --------- +# Sign flip # Find a good template subject to align other subjects to template = source_recon.find_template_subject( @@ -123,24 +120,19 @@ if False: - # ------------------------------------------------------------- - # %% Copy sf files to a single directory (makes it easier to copy minimal + # -------------------------------------------------------------------- + # Copy sf files to a single directory (makes it easier to copy minimal # files to, e.g. BMRC, for downstream analysis) os.makedirs(op.join(recon_dir, "sflip_data"), exist_ok=True) for subject, sflip_parc_file in zip(subjects, sflip_parc_files): - sflip_parc_file_from = op.join(recon_dir, subject, sflip_parc_file) - sflip_parc_file_to = op.join( recon_dir, "sflip_data", subject + "_sflip_parc-raw.fif" ) - os.system("cp -f {} {}".format(sflip_parc_file_from, sflip_parc_file_to)) - #### - if False: # Use this to copy minimum report files necessary for group workshop practical diff --git a/examples/wakeman_henson/first_level_glm.py b/examples/wakeman_henson/2023_workshop/3_first_level_glm.py similarity index 96% rename from examples/wakeman_henson/first_level_glm.py rename to examples/wakeman_henson/2023_workshop/3_first_level_glm.py index 67bd1252..1f9eaa22 100644 --- a/examples/wakeman_henson/first_level_glm.py +++ b/examples/wakeman_henson/2023_workshop/3_first_level_glm.py @@ -1,7 +1,4 @@ -#!/usr/bin/env python - -""" -Run group analysis on data from the Wakeman-Henson dataset. +"""Run group analysis on data from the Wakeman-Henson dataset. """ @@ -30,8 +27,8 @@ sessions_to_do = np.arange(0, nsessions) subj_sess_2exclude = np.zeros([nsubjects, nsessions]).astype(bool) -# ------------------------------------------------------------- -# %% Setup file names +# ---------------- +# Setup file names preproc_fif_files = [] input_fif_files = [] @@ -93,9 +90,9 @@ if not os.path.isdir(glm_subj_dir): os.makedirs(glm_subj_dir) - -# ------------------- +# ------------------------------- # Epoch first-level design matrix + for preproc_fif_file, input_fif_file, epoch_fif_file, abs_epoch_fif_file \ in zip(preproc_fif_files, input_fif_files, epoch_fif_files, abs_epoch_fif_files): @@ -117,8 +114,6 @@ epochs.load_data() epochs.save(epoch_fif_file, overwrite=True) - #### hilb - # Do hilbert transform hilb_raw = raw.copy() hilb_raw.load_data() @@ -171,8 +166,9 @@ file_to = op.join(workshop_recon_dir, subject + '/') os.system("cp -f {} {}".format(file_from, file_to)) -# ------------------- +# ------------------------------- # Setup first-level design matrix + print("\nSetting up design matrix") DC = glm.design.DesignConfig() @@ -241,8 +237,9 @@ print(DC.to_yaml()) -# ------------- +# ------------------- # Fit first-level GLM + use_hilbert = True if use_hilbert: @@ -269,9 +266,9 @@ # Create Design Matrix des = DC.design_from_datainfo(data.info) - # ------------------------------------------------------ - + # --------- # Fit Model + print("Fitting GLM") model = glm.fit.OLSModel(des, data) @@ -284,4 +281,3 @@ out.close() np.save(glm_time_file, epochs.times) - diff --git a/examples/wakeman_henson/group_level_glm.py b/examples/wakeman_henson/2023_workshop/4_group_level_glm.py similarity index 93% rename from examples/wakeman_henson/group_level_glm.py rename to examples/wakeman_henson/2023_workshop/4_group_level_glm.py index 8e3f9cdb..9a8f3a95 100644 --- a/examples/wakeman_henson/group_level_glm.py +++ b/examples/wakeman_henson/2023_workshop/4_group_level_glm.py @@ -1,7 +1,4 @@ -#!/usr/bin/env python - -""" -Run group analysis on data on the Wakeman-Henson dataset. +"""Run group analysis on data on the Wakeman-Henson dataset. """ @@ -28,8 +25,8 @@ baseline_correct = True rectify = True -# ------------------------------------------------------------- -# %% Setup file names +# ---------------- +# Setup file names glm_model_files = [] glm_time_files = [] @@ -61,8 +58,9 @@ # store which subject this belongs to subj_indices.append(sub) -# ------------------- +# -------------------------------------------- # Setup group-level design matrix in GLM tools + print("\nSetting up group design matrix") design_matrix = np.zeros([len(subj_indices), len(set(subj_indices))]) @@ -117,9 +115,9 @@ # Create GLM data data = glm.data.TrialGLMData(data=data) - # ------------------------------------------------------ - + # --------- # Fit Model + print("Fitting GLM") model = glm.fit.OLSModel(des, data) @@ -132,6 +130,3 @@ out.close() np.save(group_glm_time_file, epochs_times) - - - diff --git a/examples/wakeman_henson/plot_group_parcel_stats.py b/examples/wakeman_henson/2023_workshop/5_plot_group_parcel_stats.py similarity index 80% rename from examples/wakeman_henson/plot_group_parcel_stats.py rename to examples/wakeman_henson/2023_workshop/5_plot_group_parcel_stats.py index 05feaf2f..ad75e766 100644 --- a/examples/wakeman_henson/plot_group_parcel_stats.py +++ b/examples/wakeman_henson/2023_workshop/5_plot_group_parcel_stats.py @@ -1,7 +1,4 @@ -#!/usr/bin/env python - -""" -Run group analysis on data on the Wakeman-Henson dataset. +"""Run group analysis on data on the Wakeman-Henson dataset. """ @@ -9,13 +6,15 @@ import os import os.path as op + import numpy as np -import matplotlib.pyplot as plt -from anamnesis import obj_from_hdf5file import nibabel as nib -from osl.source_recon import parcellation, rhino +import matplotlib.pyplot as plt from nilearn import plotting +from anamnesis import obj_from_hdf5file + import glmtools as glm +from osl.source_recon import parcellation, rhino subjects_dir = "/ohba/pi/mwoolrich/datasets/WakemanHenson/ds117" subjects_dir = "/Users/woolrich/homedir/vols_data/WakeHen" @@ -31,7 +30,7 @@ show_plots = True -# ------------------------------------------------------ +# ------------------ # Load group GLM fit group_glm_model_file = op.join( @@ -45,7 +44,7 @@ model = obj_from_hdf5file(group_glm_model_file, "model") epochs_times = np.load(group_glm_time_file) -# ------------------------------------------------------ +# ------------------------------------------------------------------------------ # Output parcel-wise stats as 3D nii file in MNI space at time point of interest group_contrast = 0 @@ -62,7 +61,9 @@ print(cope_map.shape) -nii = parcellation.convert2niftii(cope_map, parcellation.find_file(parcellation_file), parcellation.find_file(mask_file)) +nii = parcellation.convert2niftii( + cope_map, parcellation.find_file(parcellation_file), parcellation.find_file(mask_file) +) cope_fname = op.join( stats_dir, @@ -102,13 +103,19 @@ ) plt.xlabel("time (s)") plt.ylabel("abs(cope)") -plt.savefig(op.join(stats_dir, 'gc{}_flc{}_tc_parc{}.png'.format(group_contrast, first_level_contrast, parcel_ind))) +plt.savefig( + op.join(stats_dir, + 'gc{}_flc{}_tc_parc{}.png'.format( + group_contrast, first_level_contrast, parcel_ind + ) + ) +) if show_plots: plt.show() -# ------------------------------------------------------------------ -# plot subject-specific copes for the first-level contrast +# -------------------------------------------------------- +# Plot subject-specific copes for the first-level contrast first_level_data = obj_from_hdf5file(group_glm_model_file, "data").data # nsess x nparcels x ntpts @@ -122,25 +129,33 @@ ) plt.xlabel("time (s)") plt.ylabel("abs(cope)") -plt.savefig(op.join(stats_dir, 'subjects_gc{}_flc{}_tc_parc{}.png'.format(group_contrast, first_level_contrast, parcel_ind))) +plt.savefig( + op.join( + stats_dir, + 'subjects_gc{}_flc{}_tc_parc{}.png'.format( + group_contrast, first_level_contrast, parcel_ind + ) + ) +) if show_plots: plt.show() -# ------------------------------------------------------------------ +# ------------------------------------------------------ # Write stats as 4D niftii file in MNI space. # Note, 4th dimension is timepoint within an epoch/trial - time_inds = np.where((epochs_times>-0.25) & (epochs_times<1))[0] cope_map = model.copes[group_contrast, :, :].T cope_map.shape -nii = parcellation.convert2niftii(cope_map, - parcellation.find_file(parcellation_file), - parcellation.find_file(mask_file), - tres=epochs_times[1]-epochs_times[0], - tmin=epochs_times[0]) +nii = parcellation.convert2niftii( + cope_map, + parcellation.find_file(parcellation_file), + parcellation.find_file(mask_file), + tres=epochs_times[1]-epochs_times[0], + tmin=epochs_times[0], +) cope_fname = op.join( stats_dir, diff --git a/examples/wakeman_henson/full_preproc/1_preprocess.py b/examples/wakeman_henson/full_preproc/1_preprocess.py new file mode 100644 index 00000000..d74cc1e6 --- /dev/null +++ b/examples/wakeman_henson/full_preproc/1_preprocess.py @@ -0,0 +1,40 @@ +"""Preprocessing. + +""" + +from dask.distributed import Client + +from osl import preprocessing, utils + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + client = Client(n_workers=6, threads_per_worker=1) + + config = """ + preproc: + - set_channel_types: {EEG061: eog, EEG062: eog, EEG063: ecg} + - filter: {l_freq: 0.5, h_freq: 125, method: iir, iir_params: {order: 5, ftype: butter}} + - notch_filter: {freqs: 50 100} + - resample: {sfreq: 250} + - bad_segments: {segment_len: 500, picks: mag} + - bad_segments: {segment_len: 500, picks: grad} + - bad_segments: {segment_len: 500, picks: mag, mode: diff} + - bad_segments: {segment_len: 500, picks: grad, mode: diff} + - ica_raw: {picks: meg, n_components: 40} + - ica_autoreject: {picks: meg, ecgmethod: correlation, eogthreshold: auto} + - interpolate_bads: {} + """ + + raw_dir = "/well/woolrich/projects/wakeman_henson/ds117" + for sub in range(1, 20): + inputs = [] + for run in range(1, 7): + inputs.append(f"{raw_dir}/sub{sub:03d}/MEG/run_{run:02d}_sss.fif") + preproc_dir = f"/well/woolrich/projects/wakeman_henson/spring23/preproc/sub{sub:02d}" + preprocessing.run_proc_batch( + config, + inputs, + outdir=preproc_dir, + overwrite=True, + dask_client=True, + ) diff --git a/examples/wakeman_henson/full_preproc/2_source_reconstruct.py b/examples/wakeman_henson/full_preproc/2_source_reconstruct.py new file mode 100644 index 00000000..a3f43964 --- /dev/null +++ b/examples/wakeman_henson/full_preproc/2_source_reconstruct.py @@ -0,0 +1,74 @@ + +"""Source reconstruction. + +This includes: coregistration, beamforming and parcellation. +""" + +import numpy as np +from dask.distributed import Client + +from osl import source_recon, utils + +source_recon.setup_fsl("/well/woolrich/projects/software/fsl") + +def fix_headshape_points(src_dir, subject, preproc_file, smri_file, epoch_file): + filenames = source_recon.rhino.get_coreg_filenames(src_dir, subject) + + # Load saved headshape and nasion files + hs = np.loadtxt(filenames["polhemus_headshape_file"]) + nas = np.loadtxt(filenames["polhemus_nasion_file"]) + lpa = np.loadtxt(filenames["polhemus_lpa_file"]) + rpa = np.loadtxt(filenames["polhemus_rpa_file"]) + + # Remove headshape points on the nose + remove = np.logical_and(hs[1] > max(lpa[1], rpa[1]), hs[2] < nas[2]) + hs = hs[:, ~remove] + + # Overwrite headshape file + utils.logger.log_or_print(f"overwritting {filenames['polhemus_headshape_file']}") + np.savetxt(filenames["polhemus_headshape_file"], hs) + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + client = Client(n_workers=16, threads_per_worker=1) + + config = """ + source_recon: + - extract_fiducials_from_fif: {} + - fix_headshape_points: {} + - compute_surfaces_coregister_and_forward_model: + include_nose: False + use_nose: False + use_headshape: True + model: Single Layer + - beamform_and_parcellate: + freq_range: [1, 45] + chantypes: [mag, grad] + rank: {meg: 60} + parcellation_file: fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz + method: spatial_basis + orthogonalisation: symmetric + """ + + raw_dir = f"/well/woolrich/projects/wakeman_henson/ds117" + src_dir = "/well/woolrich/projects/wakeman_henson/spring23/src" + + subjects = [] + preproc_files = [] + smri_files = [] + for sub in range(1, 20): + preproc_dir = f"/well/woolrich/projects/wakeman_henson/spring23/preproc/sub{sub:02d}" + for run in range(1, 7): + subjects.append(f"sub{sub:02d}_run{run:02d}") + preproc_files.append(f"{preproc_dir}/run_{run:02d}_sss/run_{run:02d}_sss_preproc_raw.fif") + smri_files.append(f"{raw_dir}/sub{sub:03d}/anatomy/highres001.nii.gz") + + source_recon.run_src_batch( + config, + src_dir=src_dir, + subjects=subjects, + preproc_files=preproc_files, + smri_files=smri_files, + extra_funcs=[fix_headshape_points], + dask_client=True, + ) diff --git a/examples/wakeman_henson/full_preproc/3_sign_flip.py b/examples/wakeman_henson/full_preproc/3_sign_flip.py new file mode 100644 index 00000000..335bef53 --- /dev/null +++ b/examples/wakeman_henson/full_preproc/3_sign_flip.py @@ -0,0 +1,46 @@ +"""Fix the dipole sign ambiguity. + +""" + +from glob import glob +from dask.distributed import Client + +from osl import source_recon, utils + +source_recon.setup_fsl("/well/woolrich/projects/software/fsl") + +# Directory containing source reconstructed data +src_dir = "/well/woolrich/projects/wakeman_henson/spring23/src" +src_files = sorted(glob(src_dir + "/*/parc/parc-raw.fif")) + +if __name__ == "__main__": + utils.logger.set_up(level="INFO") + + # Get subjects + subjects = [] + for path in src_files: + subject = path.split("/")[-3] + subjects.append(subject) + + # Find a good template subject to match others to + template = source_recon.find_template_subject( + src_dir, subjects, n_embeddings=15, standardize=True + ) + + # Settings + config = f""" + source_recon: + - fix_sign_ambiguity: + template: {template} + n_embeddings: 15 + standardize: True + n_init: 3 + n_iter: 3000 + max_flips: 20 + """ + + # Setup parallel processing + client = Client(n_workers=16, threads_per_worker=1) + + # Run sign flipping + source_recon.run_src_batch(config, src_dir, subjects, dask_client=True) diff --git a/examples/wakeman_henson/full_preproc/4_save_npy.py b/examples/wakeman_henson/full_preproc/4_save_npy.py new file mode 100644 index 00000000..7515a4e2 --- /dev/null +++ b/examples/wakeman_henson/full_preproc/4_save_npy.py @@ -0,0 +1,11 @@ +"""Save training data as npy files. + +""" + +from glob import glob +from osl_dynamics.data import Data + +files = sorted(glob("/well/woolrich/projects/wakeman_henson/spring23/src/*/sflip_parc-raw.fif")) +data = Data(files, picks="misc", reject_by_annotation="omit") +data.save("/well/woolrich/projects/wakeman_henson/spring23/training") +data.delete_dir() diff --git a/examples/wakeman_henson/wakeman_henson.py b/examples/wakeman_henson/run_wakeman_henson.py similarity index 97% rename from examples/wakeman_henson/wakeman_henson.py rename to examples/wakeman_henson/run_wakeman_henson.py index f6bd5bd6..a9453c56 100644 --- a/examples/wakeman_henson/wakeman_henson.py +++ b/examples/wakeman_henson/run_wakeman_henson.py @@ -1,9 +1,5 @@ -#!/usr/bin/env python - -"""Run RHINO-based source reconstruction on the Wakeman-Henson dataset. - -This script can be copy and pasted into any directory on hbaws and executed with: -python wakeman_henson.py +"""Run RHINO-based source reconstruction on the Wakeman-Henson dataset and do +first-level GLM analysis. See run_wakeman_henson.m for matlab equivalent. """ @@ -129,6 +125,7 @@ # ------------------- # Setup design matrix + print("\nSetting up design matrix") DC = glm.design.DesignConfig() @@ -216,9 +213,9 @@ epochs.pick(chantypes) data = glm.io.load_mne_epochs(epochs) - # ------------------------------------------------------ - + # --------- # Fit Model + print("Fitting GLM") model_sensor = glm.fit.OLSModel(des, data) @@ -231,9 +228,9 @@ model_sensor.to_hdf5(out.create_group("model_sensor")) out.close() - # ------------------------------------------------------ - + # ---------------- # Load subject GLM + model_sensor = obj_from_hdf5file(glmname, "model_sensor") # Make MNE object with contrast @@ -252,6 +249,7 @@ # ------------------------------------------------------------------ # Compute info for source reconstruction (coreg, forward model, BEM) + print("Computing info for source reconstruction") subjects_dir = op.join(out_dir, "coreg") @@ -332,6 +330,7 @@ # ------------------------------------------- # Apply source reconstruction to epoched data + print("Applying source reconstruction to epoched data") # Epoch the preprocessed data again from scratch and drop bad epochs @@ -375,6 +374,7 @@ # ------------------------------------ # Fit GLM to source reconstructed data + print("Source Space Analysis") # Turns this into a ntrials x nsources x ntpts array @@ -390,9 +390,9 @@ print("Please close the pop-up figures to continue running the code") des.plot_summary(show=True, savepath=preproc_fif_file.replace(".fif", "_design.png")) -# ------------------------------------------------------ - +# --------- # Fit Model + print("Fitting GLM") model_source = glm.fit.OLSModel(des, data) @@ -407,6 +407,7 @@ # -------------------------------------- # Compute stats of interest from GLM fit + print("Computing stats from GLM fit") acopes = [] @@ -429,6 +430,7 @@ # ------------------------------------------------------ # Output stats as 3D nii files at time point of interest + print("Saving stats nii") stats_dir = op.join(subjects_dir, subject, "rhino", "stats") @@ -477,6 +479,7 @@ # -------------------------------------------------------------------------- # Convert cope to standard brain grid in MNI space, for doing group stats # (sourcespace_epoched_data and therefore acopes are in head/polhemus space) + print("Converting COPE to MNI space") ( diff --git a/osl/source_recon/beamforming.py b/osl/source_recon/beamforming.py index aa6ddea1..8800922a 100644 --- a/osl/source_recon/beamforming.py +++ b/osl/source_recon/beamforming.py @@ -43,6 +43,7 @@ from osl.source_recon.rhino import utils as rhino_utils from osl.utils.logger import log_or_print + def get_beamforming_filenames(subjects_dir, subject): """Get beamforming filenames. @@ -644,7 +645,18 @@ def _make_lcmv( # Compute spatial filter n_orient = 3 if is_free_ori else 1 W, max_power_ori = _compute_beamformer( - G, Cm, reg, n_orient, weight_norm, pick_ori, reduce_rank, rank_int, inversion=inversion, nn=nn, orient_std=orient_std, whitener=whitener, + G, + Cm, + reg, + n_orient, + weight_norm, + pick_ori, + reduce_rank, + rank_int, + inversion=inversion, + nn=nn, + orient_std=orient_std, + whitener=whitener, ) # Get src type to store with filters for _make_stc @@ -961,13 +973,13 @@ def _prepare_beamformer_input( Modified version of mne.beamformer._prepare_beamformer_input. """ - # Lines marked MWW or CG for where code has been changed. + # Lines marked MWW are where code has been changed. # MWW # _check_option('pick_ori', pick_ori, ('normal', 'max-power', 'vector', None)) _check_option("pick_ori", pick_ori, ("normal", "max-power", "vector", "max-power-pre-weight-norm", None)) - # MWW, CG + # MWW # Restrict forward solution to selected vertices #if label is not None: # _, src_sel = label_src_vertno_sel(label, forward["src"]) @@ -976,7 +988,7 @@ def _prepare_beamformer_input( if loose is None: loose = 0.0 if is_fixed_orient(forward) else 1.0 - # MWW, CG + # MWW #if noise_cov is None: # noise_cov = make_ad_hoc_cov(info, std=1.0) @@ -1087,11 +1099,8 @@ def voxel_timeseries( log_or_print(f"bandpass filtering: {freq_range}") data = data.filter(l_freq=freq_range[0], h_freq=freq_range[1], method="iir", iir_params={"order": 5, "ftype": "butter"}) - # Load the beamformer + # Load the beamformer and apply filters = load_lcmv(subjects_dir, subject) - - # Apply the beamformer - log_or_print("applying beamformer") bf_data = apply_lcmv(data, filters, reject_by_annotation) # Transform to MNI space diff --git a/osl/source_recon/parcellation/parcellation.py b/osl/source_recon/parcellation/parcellation.py index 11fb7ed4..eaed6a9f 100644 --- a/osl/source_recon/parcellation/parcellation.py +++ b/osl/source_recon/parcellation/parcellation.py @@ -68,11 +68,13 @@ def parcellate_timeseries(parcellation_file, voxel_timeseries, voxel_coords, met Parameters ---------- + parcellation_file : str + Parcellation file (or path to parcellation file). voxel_timeseries : numpy.ndarray - nvoxels x ntpts, or nvoxels x ntpts x ntrials Data to be parcellated. Data is assumed to be in same space as the parcellation (e.g. typically corresponds - to the output from rhino.resample_recon_ts). + (nvoxels x ntpts) or (nvoxels x ntpts x ntrials) data to be parcellated. + Data is assumed to be in same space as the parcellation (e.g. typically corresponds to the output from beamforming.transform_recon_timeseries). voxel_coords : numpy.ndarray - (nvoxels x 3) coordinates of voxel_timeseries in mm in same space as parcellation (e.g. typically corresponds to the output from rhino.resample_recon_ts). + (nvoxels x 3) coordinates of voxel_timeseries in mm in same space as parcellation (e.g. typically corresponds to the output from beamforming.transform_recon_timeseries). method : str 'pca' - take 1st PC in each parcel 'spatial_basis' - The parcel time-course for each spatial map is the 1st PC from all voxels, weighted by the spatial map. @@ -98,10 +100,10 @@ def _get_parcel_timeseries(voxel_timeseries, parcellation_asmatrix, method="spat Parameters ---------- - parcellation_asmatrix: numpy.ndarray - nvoxels x nparcels and is assumed to be on the same grid as voxel_timeseries voxel_timeseries : numpy.ndarray - nvoxels x ntpts or nvoxels x ntpts x ntrials and is assumed to be on the same grid as parcellation (typically output by rhino.resample_recon_ts) + (nvoxels x ntpts) or (nvoxels x ntpts x ntrials) and is assumed to be on the same grid as parcellation (typically output by beamforming.transform_recon_timeseries). + parcellation_asmatrix: numpy.ndarray + (nvoxels x nparcels) and is assumed to be on the same grid as voxel_timeseries. method : str 'pca' - take 1st PC of voxels 'spatial_basis' - The parcel time-course for each spatial map is the 1st PC from all voxels, weighted by the spatial map. @@ -269,10 +271,10 @@ def _get_parcel_timeseries(voxel_timeseries, parcellation_asmatrix, method="spat def _resample_parcellation(parcellation_file, voxel_coords, working_dir=None): - """Resample parcellation so that its voxel coords correspond (using nearest neighbour) to passed in voxel_coords. Passed in voxel_coords - and parcellation must be in the same space, e.g. MNI. + """Resample parcellation so that its voxel coords correspond (using nearest neighbour) to passed in voxel_coords. + Passed in voxel_coords and parcellation must be in the same space, e.g. MNI. - Used to make sure that the parcellation's voxel coords are the same as the voxel coords for some timeseries data, before calling get_parcel_timeseries. + Used to make sure that the parcellation's voxel coords are the same as the voxel coords for some timeseries data, before calling _get_parcel_timeseries. Parameters ---------- @@ -292,16 +294,16 @@ def _resample_parcellation(parcellation_file, voxel_coords, working_dir=None): log_or_print(f"gridstep = {gridstep} mm") parcellation_file = find_file(parcellation_file) - pth, parcellation_name = op.split(op.splitext(op.splitext(parcellation_file)[0])[0]) + path, parcellation_name = op.split(op.splitext(op.splitext(parcellation_file)[0])[0]) if working_dir is None: - working_dir = pth + working_dir = path else: os.makedirs(working_dir, exist_ok=True) parcellation_resampled = op.join(working_dir, parcellation_name + "_{}mm.nii.gz".format(gridstep)) - # Create std brain of the required resolution + # Create standard brain of the required resolution # # Command: flirt -in -ref -out -applyisoxfm # @@ -566,17 +568,17 @@ def _parcel_timeseries2nii( parcel_timeseries_data: numpy.ndarray Needs to be nparcels x ntpts voxel_weightings : numpy.ndarray - nvoxels x nparcels voxel weightings for each parcel to compute parcel_timeseries from voxel_timeseries. + (nvoxels x nparcels) voxel weightings for each parcel to compute parcel_timeseries from voxel_timeseries. voxel_assignments : bool numpy.ndarray - nvoxels x nparcels, Boolean assignments indicating for each voxel the winner takes all parcel it belongs to. + (nvoxels x nparcels) boolean assignments indicating for each voxel the winner takes all parcel it belongs to. voxel_coords : numpy.ndarray - (nvoxels x 3) coordinates of voxel_timeseries in mm in same space as parcellation (e.g. typically corresponds to the output from rhino.resample_recon_ts). + (nvoxels x 3) coordinates of voxel_timeseries in mm in same space as parcellation (e.g. typically corresponds to the output from beamforming.transform_recon_timeseries). working_dir : str Directory name to put files in. out_nii_fname : str Output name to put files in. - times : (ntpts, ) numpy.ndarray - Times points in seconds. Will assume that time points are regularly spaced. Used to set nii file up correctly. + times : array + (ntpts,) times points in seconds. Will assume that time points are regularly spaced. Used to set nii file up correctly. method : str "weights" or "assignments" @@ -586,10 +588,10 @@ def _parcel_timeseries2nii( Output nii filename, will be output at spatial resolution of parcel_timeseries['voxel_coords']. """ parcellation_file = find_file(parcellation_file) - pth, parcellation_name = op.split(op.splitext(op.splitext(parcellation_file)[0])[0]) + path, parcellation_name = op.split(op.splitext(op.splitext(parcellation_file)[0])[0]) if working_dir is None: - working_dir = pth + working_dir = path if out_nii_fname is None: out_nii_fname = op.join(working_dir, parcellation_name + "_timeseries.nii.gz") @@ -615,7 +617,7 @@ def _parcel_timeseries2nii( gridstep = int(rhino_utils.get_gridstep(voxel_coords.T) / 1000) # Sample parcellation_mask to the desired resolution - pth, ref_brain_name = op.split(op.splitext(op.splitext(parcellation_mask_file)[0])[0]) + path, ref_brain_name = op.split(op.splitext(op.splitext(parcellation_mask_file)[0])[0]) parcellation_mask_resampled = op.join(working_dir, ref_brain_name + "_{}mm_brain.nii.gz".format(gridstep)) # create std brain of the required resolution @@ -659,7 +661,7 @@ def convert2niftii(parc_data, parcellation_file, mask_file, tres=1, tmin=0): parc_data : np.ndarray (nparcels) or (nvolumes x nparcels) parcel data. parcellation_file : str - Path to niftii parcellation file.. + Path to niftii parcellation file. mask_file : str Path to niftii parcellation mask file. tres : float @@ -736,7 +738,7 @@ def convert2mne_raw(parc_data, raw, parcel_names=None, extra_chans="stim"): Parameters ---------- parc_data : np.ndarray - nparcels x ntpts parcel data + (nparcels x ntpts) parcel data. raw : mne.Raw mne.io.raw object that produced parc_data via source recon and parcellation. Info such as timings and bad segments will be copied from this to parc_raw. parcel_names : list of str @@ -802,10 +804,10 @@ def convert2mne_epochs(parc_data, epochs, parcel_names=None): Parameters ---------- parc_data : np.ndarray - nparcels x ntpts x epochs parcel data + (nparcels x ntpts x epochs) parcel data. epochs : mne.Epochs mne.io.raw object that produced parc_data via source recon and parcellation. Info such as timings and bad segments will be copied from this to parc_raw. - parcel_names : list(str) + parcel_names : list of str List of strings indicating names of parcels. If None then names are set to be parcel_0,...,parcel_{n_parcels-1}. Returns diff --git a/osl/source_recon/rhino/utils.py b/osl/source_recon/rhino/utils.py index d2ca429c..fd88e695 100644 --- a/osl/source_recon/rhino/utils.py +++ b/osl/source_recon/rhino/utils.py @@ -1240,7 +1240,7 @@ def extract_rhino_files(old_subjects_dir, new_subjects_dir, subjects="all", excl else: raise FileNotFoundError(old_file) - # Special case + # Special case std_brains = glob(f"{old_dir}/MNI152_T1_*_brain.nii.gz") for std_brain in std_brains: copy(std_brain, std_brain.replace(old_dir, new_dir))