Merge pull request #356 from OHBA-analysis/osl-toolbox-paper

added toolbox paper scripts

examples/toolbox-paper/1_preprocessing.py

@@ -0,0 +1,56 @@
+import os
+from dask.distributed import Client
+
+import osl_ephys
+
+
+if __name__ == "__main__": 
+  client = Client(n_workers=16, threads_per_worker=1) # specify to enable parallel processing
+  basedir = "ds117"
+
+  config = """
+    meta:
+      event_codes:
+        famous/first: 5
+        famous/immediate: 6
+        famous/last: 7
+        unfamiliar/first: 13
+        unfamiliar/immediate: 14
+        unfamiliar/last: 15
+        scrambled/first: 17
+        scrambled/immediate: 18
+        scrambled/last: 19
+    preproc:
+      - find_events: {min_duration: 0.005}
+      - 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}
+      - bad_channels: {picks: mag, significance_level: 0.1}
+      - bad_channels: {picks: grad, significance_level: 0.1}
+      - ica_raw: {picks: meg, n_components: 40}
+      - ica_autoreject: {picks: meg, ecgmethod: correlation, eogmethod: correlation,
+            eogthreshold: 0.35, apply: False}
+      - interpolate_bads: {reset_bads: False}
+      """
+
+  # Study utils enables selection of existing paths using various wild cards
+  study = osl_ephys.utils.Study(os.path.join(basedir, "sub{sub_id}/MEG/run_{run_id}_raw.fif"))
+  inputs = sorted(study.get())
+
+  # specify session names and output directory
+  subjects = [f"sub{i+1:03d}-run{j+1:02d}" for i in range(19) for j in range(6)]
+  outdir = os.path.join(basedir, "processed")
+
+  osl_ephys.preprocessing.run_proc_batch(
+    config,
+    inputs,
+    subjects,
+    outdir,
+    dask_client=True,
+    random_seed=2280431064,
+  )

examples/toolbox-paper/2_source-reconstruct.py

@@ -0,0 +1,68 @@
+import os
+import numpy as np
+from dask.distributed import Client
+from osl_ephys import source_recon, utils
+
+source_recon.setup_fsl("~/fsl") # FSL needs to be installed
+
+def fix_headshape_points(outdir, subject):
+    filenames = source_recon.rhino.get_coreg_filenames(outdir, 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_polhemus_from_info: {}
+        - fix_headshape_points: {}
+        - compute_surfaces:
+            include_nose: False
+        - coregister:
+            use_nose: False
+            use_headshape: True
+        - forward_model:
+            model: Single Layer
+        - beamform_and_parcellate:
+            freq_range: [1, 45]
+            chantypes: [mag, grad]
+            rank: {meg: 60}
+            parcellation_file: Glasser52_binary_space-MNI152NLin6_res-8x8x8.nii.gz
+            method: spatial_basis
+            orthogonalisation: symmetric
+    """
+
+    basedir = "ds117"
+    proc_dir = os.path.join(basedir, "processed")
+
+    # Define inputs
+    subjects = [f"sub{i+1:03d}-run{j+1:02d}" for i in range(19) for j in range(6)]
+    preproc_files = sorted(utils.Study(os.path.join(proc_dir, "sub{sub_id}-    run{run_id}/sub{sub_id}-run{run_id}_preproc-raw.fif")).get())
+    smri_files = np.concatenate([[smri_file]*6 for smri_file in sorted(utils.Study(os.path.join(basedir, "sub{sub_id}/anatomy/highres001.nii.gz"))).get()])
+
+    # Run source batch
+    source_recon.run_src_batch(
+        config,
+        outdir=proc_dir,
+        subjects=subjects,
+        preproc_files=preproc_files,
+        smri_files=smri_files,
+        extra_funcs=[fix_headshape_points],
+        dask_client=True,
+        random_seed=1392754308,
+    )
+

examples/toolbox-paper/3_sign-flip.py

@@ -0,0 +1,46 @@
+
+import os
+from glob import glob
+from dask.distributed import Client
+
+from osl_ephys import source_recon, utils
+
+source_recon.setup_fsl("~/fsl")
+
+# Directory containing source reconstructed data
+proc_dir = "ds117/processed"
+src_files = sorted(utils.Study(os.path.join(proc_dir, 
+"sub{sub_id}-run{run_id}/parc/parc-raw.fif")).get())
+
+if __name__ == "__main__":
+    utils.logger.set_up(level="INFO")
+
+    subjects = [f"sub{i+1:03d}-run{j+1:02d}" for i in range(19) for j in range(6)]
+
+    # Find a good template subject to match others to
+    template = source_recon.find_template_subject(
+        proc_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, proc_dir, subjects, dask_client=True, random_seed=3116145039)
+
+
+
+
+

examples/toolbox-paper/4_stats.py

@@ -0,0 +1,46 @@
+import os
+import numpy as np
+import glmtools
+import matplotlib.pyplot as plt
+from dask.distributed import Client
+from osl_ephys import preprocessing, glm
+
+
+if __name__ == "__main__":
+    client = Client(n_workers=16, threads_per_worker=1)
+
+    config = """
+      preproc:
+        - read_dataset: {ftype: sflip_parc-raw}
+        - epochs: {picks: misc, tmin: -0.2, tmax: 0.3}
+        - glm_add_regressor: {name: famous, rtype: Categorical, codes: [5 6 7]}
+        - glm_add_regressor: {name: unfamiliar, rtype: Categorical, codes: [13 14 15]}
+        - glm_add_regressor: {name: scrambled, rtype: Categorical, codes: [17 18 19]}
+        - glm_add_contrast: {name: Mean, values: {famous: 1/3, unfamiliar: 1/3, scrambled: 1/3}}
+        - glm_add_contrast: {name: Faces-Scrambled, values: {famous: 1, unfamiliar: 1, scrambled: -2}}
+        - glm_fit: {target: epochs, method: glm_epochs}
+      group:
+        - glm_add_regressor: {name: Subject, rtype: Categorical, key: Subject, codes: unique}
+        - glm_add_contrast: {name: Mean, values: unique, key: Subject}
+        - glm_fit: {method: epochs, tmin: 0.05, tmax: 0.3}
+        - glm_permutations: {method: epochs, target: group_glm, contrast: Mean, type: max, nperms: 1000, threshold: 0.99}
+    """    
+    proc_dir = "ds117/processed"
+    src_files = sorted(utils.Study(os.path.join(proc_dir, 
+"sub{sub_id}-run{run_id}", "sub{sub_id}-run{run_id}_sflip_parc-raw.fif")).get())    
+    subjects = [f"sub{i+1:03d}-run{j+1:02d}" for i in range(19) for j in range(6)]
+    covs = [f"Subject": [sub.split("-")[0]for sub in subjects]
+
+    preprocessing.run_proc_batch(
+        config,
+        src_files,
+        subjects,
+        outdir=proc_dir,
+        ftype='raw',
+        covs=covs,
+        dask_client=True,
+        overwrite=True,
+        gen_report=False,
+        skip_save=['events', 'raw', 'ica', 'event_id', 'sflip_parc-raw'],
+    )
+

examples/toolbox-paper/README.md

@@ -0,0 +1,22 @@
+# OSL Toolbox Paper
+
+The scripts for the toolbox paper can be found here in the order they appear in the manuscript. You can download the data from [OpenfMRI](https://openfmri.org/s3-browser/?prefix=ds000117/ds000117_R0.1.1/compressed/). Extract the `tar.gz` files in a folder called `ds117`. Note: you may have to change the base directory in the scripts to match the directory where you store the data. You also need to [install FSL](https://fsl.fmrib.ox.ac.uk/fsl/docs/#/install/index), and make sure that `source_recon.setup_fsl()` in `2_source-reconstruct.py` and `3_sign-flip.py` is pointing to the correct directory.
+
+If you wish to fully reproduce the analysis pipeline install the environment specified in `osl-toolbox-paper.yml`, and set `random_seed` in `run_proc_batch` and `run_src_batch` according to the seed found in the logfiles on [OSF](https://osf.io/2rnyg/).
+
+## Manual preprocessing
+
+Note that in the toolbox paper, automatically labeled ICA components were manually refined for the following sessions:
+- sub008-ses03
+- sub019-ses01
+- sub019-ses02
+- sub019-ses03
+- sub019-ses04
+- sub019-ses05
+- sub019-ses06
+- sub010-ses05
+
+This was done by running the following command line function iteratively, replacing "session". 
+`osl_ica_label None processed "session"`
+After all specified sessions were refined, all automatically/manually labeled components were removed from the preprocessed MEG data using the command line call
+`osl_ica_apply processed`