From 81726e4c515bce4cae3bdebbae8ce183e51271b4 Mon Sep 17 00:00:00 2001 From: Urszula Date: Wed, 5 Jul 2023 23:47:27 -0500 Subject: [PATCH 1/4] The most basic (but enough) version of QC --- coglib/meeg/qc/P00_bids_conversion.py | 7 - coglib/meeg/qc/P00_run_qc.py | 42 ---- coglib/meeg/qc/P00_run_qc_epochs.py | 82 ------- coglib/meeg/qc/QC_epochs.py | 31 --- coglib/meeg/qc/QC_processing.py | 35 +-- coglib/meeg/qc/QC_processing_eeg.py | 329 -------------------------- coglib/meeg/qc/srun_qc.py | 29 +++ 7 files changed, 31 insertions(+), 524 deletions(-) delete mode 100644 coglib/meeg/qc/P00_run_qc_epochs.py delete mode 100644 coglib/meeg/qc/QC_processing_eeg.py create mode 100644 coglib/meeg/qc/srun_qc.py diff --git a/coglib/meeg/qc/P00_bids_conversion.py b/coglib/meeg/qc/P00_bids_conversion.py index a8d8dc1..a2d32d2 100644 --- a/coglib/meeg/qc/P00_bids_conversion.py +++ b/coglib/meeg/qc/P00_bids_conversion.py @@ -6,13 +6,6 @@ https://mne.tools/mne-bids/stable/index.html -Questions/Issues: - - what to write in the participants and dataset_description metadata - files? - - what session ID should we give to the anat scan (e.g., v0, v2, mri)? - - for visit 2, what to count the replay runs? Continue the count from - where it was left from the VG (run 4) or restart from run 1? - Notes: - the conversion must be done after reading the events. Here, the event list includes all the triggers/events diff --git a/coglib/meeg/qc/P00_run_qc.py b/coglib/meeg/qc/P00_run_qc.py index ffd8939..9950518 100644 --- a/coglib/meeg/qc/P00_run_qc.py +++ b/coglib/meeg/qc/P00_run_qc.py @@ -37,46 +37,4 @@ has_eeg = True -# # ============================================================================= -# # DEFINE PREPROCESSING STEPS -# # ============================================================================= - -# def pre_step1(): -# P01_maxwell_filtering.run_maxwell_filter(subject_id, -# visit_id) -# if has_eeg: -# P02_find_bad_eeg.find_bad_eeg(subject_id, -# visit_id, -# has_eeg) -# P03_artifact_annotation.artifact_annotation(subject_id, -# visit_id, -# has_eeg, -# # threshold_muscle, -# ) -# P04_extract_events.run_events(subject_id, -# visit_id) -# P05_run_ica.run_ica(subject_id, -# visit_id, -# has_eeg) - -# def pre_step2( -# # meg_ica_eog=opt.mICA_eog, meg_ica_ecg=opt.mICA_ecg, -# # eeg_ica_eog=opt.eICA_eog, eeg_ica_ecg=opt.eICA_ecg, -# ): -# P06_apply_ica.apply_ica(subject_id, -# visit_id, -# has_eeg) - -# P07_make_epochs.run_epochs(subject_id, -# visit_id, -# has_eeg) - - -# ============================================================================= -# RUN -# ============================================================================= -# if opt.step == '1': -# pre_step1() -# elif opt.step == '2': -# pre_step2() QC_processing.run_qc_processing(subject_id, visit_id, has_eeg) diff --git a/coglib/meeg/qc/P00_run_qc_epochs.py b/coglib/meeg/qc/P00_run_qc_epochs.py deleted file mode 100644 index 703f55a..0000000 --- a/coglib/meeg/qc/P00_run_qc_epochs.py +++ /dev/null @@ -1,82 +0,0 @@ -# -*- coding: utf-8 -*- -""" -@author: Urszula Górska gorska@wisc.edu -""" - -import argparse - -import QC_epochs - -# ============================================================================= -# PARSER SETTINGS -# ============================================================================= - -parser=argparse.ArgumentParser() -parser.add_argument('--sub', type=str, default='SA101', help='subject_id') -parser.add_argument('--visit', type=str, default='V1', help='visit_id') - -opt=parser.parse_args() - -# ============================================================================= -# SESSION-SPECIFIC SETTINGS -# ============================================================================= - -subject_id = opt.sub -visit_id = opt.visit - -# Find out whether the participant has EEG data -if visit_id.upper() == 'V1': - if subject_id.upper() in ['SA101', 'SA102', 'SA103', 'SA104']: - has_eeg = False - else: - has_eeg = True -elif visit_id.upper() == 'V2': - if subject_id.upper() in ['SA104', 'SA106']: - has_eeg = False - else: - has_eeg = True - - -# # ============================================================================= -# # DEFINE PREPROCESSING STEPS -# # ============================================================================= - -# def pre_step1(): -# P01_maxwell_filtering.run_maxwell_filter(subject_id, -# visit_id) -# if has_eeg: -# P02_find_bad_eeg.find_bad_eeg(subject_id, -# visit_id, -# has_eeg) -# P03_artifact_annotation.artifact_annotation(subject_id, -# visit_id, -# has_eeg, -# # threshold_muscle, -# ) -# P04_extract_events.run_events(subject_id, -# visit_id) -# P05_run_ica.run_ica(subject_id, -# visit_id, -# has_eeg) - -# def pre_step2( -# # meg_ica_eog=opt.mICA_eog, meg_ica_ecg=opt.mICA_ecg, -# # eeg_ica_eog=opt.eICA_eog, eeg_ica_ecg=opt.eICA_ecg, -# ): -# P06_apply_ica.apply_ica(subject_id, -# visit_id, -# has_eeg) - -# P07_make_epochs.run_epochs(subject_id, -# visit_id, -# has_eeg) - - -# ============================================================================= -# RUN -# ============================================================================= -# if opt.step == '1': -# pre_step1() -# elif opt.step == '2': -# pre_step2() -QC_epochs.run_qc_epochs(subject_id, visit_id, has_eeg) diff --git a/coglib/meeg/qc/QC_epochs.py b/coglib/meeg/qc/QC_epochs.py index 31e1785..5ca3510 100644 --- a/coglib/meeg/qc/QC_epochs.py +++ b/coglib/meeg/qc/QC_epochs.py @@ -189,37 +189,6 @@ def run_qc_epochs(subject_id, visit_id, has_eeg): epochs_irr_O = epochs['Trial_type == "Filler" and Stimuli_type == "Blank"'] print(epochs_irr_O) - #types0 = ['Filler', 'Probe'] - #type1 = ['Face', 'Object', 'Blank'] - #location = ['Upper Left', 'Upper Right', 'Lower Right', 'Lower Left'] - #response = ['Seen', 'Unseen'] - - #print(epochs[['face1', 'face2','face3','face4','face5','face6','face7','face8','face9','face10']]) - #print(epochs[['object1', 'object2','object3','object4','object5','object6','object7','object8','object9','object10']]) - #print(epochs[['letter1', 'letter2','letter3','letter4','letter5','letter6','letter7','letter8','letter9','letter10']]) - #print(epochs[['false1', 'false2','false3','false4','false5','false6','false7','false8','false9','false10']]) - - # Get rejection thresholds - MEG only - #reject = get_rejection_threshold(epochs, ch_types=['mag', 'grad'], #'eeg'], - # decim=2) - - # Drop bad epochs based on peak-to-peak magnitude - #epochs.drop_bad(reject=reject) - - #print("VERY_IMPORTANT EPOCHS faces after drop") - #epochs_rel_F = epochs['Task_relevance == "Relevant non-target" and Category == "face"'] - #print(epochs_rel_F) - #print(epochs[['face1', 'face2','face3','face4','face5','face6','face7','face8','face9','face10']]) - #print(epochs[['object1', 'object2','object3','object4','object5','object6','object7','object8','object9','object10']]) - #print(epochs[['letter1', 'letter2','letter3','letter4','letter5','letter6','letter7','letter8','letter9','letter10']]) - #print(epochs[['false1', 'false2','false3','false4','false5','false6','false7','false8','false9','false10']]) - - #print("VERY_IMPORTANT DROP LOG") - #print(epochs.drop_log) - - # Plot percentage of rejected epochs per channel - #fig1 = epochs.plot_drop_log() - #pdf.savefig(fig1) plt.close() diff --git a/coglib/meeg/qc/QC_processing.py b/coglib/meeg/qc/QC_processing.py index 8e2da55..ae6a1d8 100644 --- a/coglib/meeg/qc/QC_processing.py +++ b/coglib/meeg/qc/QC_processing.py @@ -344,7 +344,7 @@ def run_qc_processing(subject_id, visit_id, has_eeg): # check=False) #epochs = mne.read_epochs(bids_path_epo.fpath, preload=False) - print("VERY_IMPORTANT EPOCHS faces") # save outputs to pdf for everyone to access + print("EPOCHS before drop check - faces") # save outputs to pdf for everyone to access epochs_rel_F = epochs['Task_relevance == "Relevant non-target" and Category == "face"'] print(epochs_rel_F) #print(epochs[['face1', 'face2','face3','face4','face5','face6','face7','face8','face9','face10']]) @@ -356,7 +356,7 @@ def run_qc_processing(subject_id, visit_id, has_eeg): # Drop bad epochs based on peak-to-peak magnitude epochs.drop_bad(reject=reject) - print("VERY_IMPORTANT AFTER DROP :)") + print("EPOCHS AFTER DROP :)") print("FACES task relevant") epochs_rel_F = epochs['Task_relevance == "Relevant non-target" and Category == "face"'] print(epochs_rel_F) @@ -423,38 +423,7 @@ def run_qc_processing(subject_id, visit_id, has_eeg): #fig1 = epochs.plot_drop_log() #pdf.savefig(fig1) #plt.close() - - - # Epoch raw data 3/3 - #epochs = mne.Epochs(raw, - # events, - # events_id, - # tmin, tmax, - # baseline=None, - # proj=True, - # picks=picks, - # detrend=1, - # reject=None, - # reject_by_annotation=True, - # verbose=True) - - #del raw - - # Add metadata - #epochs.metadata = metadata - # Get rejection thresholds - all - #reject = get_rejection_threshold(epochs, ch_types=['mag', 'grad', 'eeg'], - # decim=2) - - # Drop bad epochs based on peak-to-peak magnitude - #epochs.drop_bad(reject=reject) - - # Plot percentage of rejected epochs per channel - #fig1 = epochs.plot_drop_log() - #pdf.savefig(fig1) - #plt.close() - if __name__ == '__main__': subject_id = input("Type the subject ID (e.g., SA101)\n>>> ") visit_id = input("Type the visit ID (V1 or V2)\n>>> ") diff --git a/coglib/meeg/qc/QC_processing_eeg.py b/coglib/meeg/qc/QC_processing_eeg.py deleted file mode 100644 index 70c89a0..0000000 --- a/coglib/meeg/qc/QC_processing_eeg.py +++ /dev/null @@ -1,329 +0,0 @@ -import os -import os.path as op -import matplotlib.pyplot as plt -import pandas as pd - - -import mne -import mne_bids -from pyprep.prep_pipeline import PrepPipeline -from mne.preprocessing import annotate_muscle_zscore -from autoreject import get_rejection_threshold -from numpy import arange - -from config import (bids_root, tmin, tmax) - -from matplotlib.backends.backend_pdf import PdfPages - -from qc.maxwell_filtering import run_maxwell_filter -from qc.extract_events import run_events - -from qc.viz_psd import viz_psd - -run_part_1 = True -run_part_2 = True -run_part_3 = False - -def run_qc_processing(subject_id, visit_id, has_eeg): - # Set path to qc derivatives and create the related folders - prep_deriv_root = op.join(bids_root, "derivatives", "qc", visit_id) - if not op.exists(prep_deriv_root): - os.makedirs(prep_deriv_root) - - print("Processing subject: %s" % subject_id) - - raw_list = list() - events_list = list() - metadata_list = list() - - with PdfPages(op.join(prep_deriv_root, subject_id + '_' + visit_id + '_MEG_QC.pdf')) as pdf: - #region first page - - FirstPage = plt.figure(figsize=(8,1), dpi=108) - FirstPage.clf() - plt.axis('off') - plt.text(0.5, 0.5, subject_id, transform=FirstPage.transFigure, size=16, ha="center") - pdf.savefig(FirstPage) - plt.close() - - #endregion first page - - bids_eo_path = mne_bids.BIDSPath(root=bids_root, - datatype='meg', - subject=subject_id, - session=visit_id, - task='rest', - extension='.fif' - ) - # Read raw data - raw_eo = mne_bids.read_raw_bids(bids_eo_path) - - fig_eo = plt.figure(figsize=(9,6), dpi=108) - ax1 = plt.subplot2grid((2,2), (0,0), colspan=2) - ax1.set_title('Resting EO spectra') - ax2 = plt.subplot2grid((2,2), (1,0), colspan=2) - - raw_eo.plot_psd(fmax=100, ax=[ax1, ax2], picks = ['meg']) - pdf.savefig(fig_eo) - plt.close() - - data_path = os.path.join(bids_root, f"sub-{subject_id}", f"ses-{visit_id}", "meg") - run = 0 - for fname in os.listdir(data_path): - if fname.endswith(".fif") and "run" in fname: - run = run + 1 - print(" Run: %s" % run) - - # Set task - if 'dur' in fname: - bids_task = 'dur' - elif 'vg' in fname: - bids_task = 'vg' - elif 'replay' in fname: - bids_task = 'replay' - else: - raise ValueError("Error: could not find the task for %s" % fname) - - # Set BIDS path - bids_path = mne_bids.BIDSPath( - root=bids_root, - subject=subject_id, - datatype='meg', - task=bids_task, - run=f"{run:02}", - session=visit_id, - extension='.fif') - - # Read raw data - raw = mne_bids.read_raw_bids(bids_path) - - #region PART 2a - MEG data filtering using Maxwell filters, method - defined in config - - if run_part_2: - # Find initial head position - if run == 1: - destination = raw.info['dev_head_t'] - - # # Set BIDS path - # bids_path_sss = mne_bids.BIDSPath( - # root=op.join(bids_root, "derivatives", "preprocessing"), - # subject=subject_id, - # datatype='meg', - # task=bids_task, - # run=f"{run:02}", - # session=visit_id, - # suffix="sss", - # extension='.fif', - # check=False) - - # # Read raw data - # raw_sss = mne_bids.read_raw_bids(bids_path_sss).load_data() - - raw_sss, bad_chan = run_maxwell_filter(raw, destination, bids_path.meg_crosstalk_fpath, bids_path.meg_calibration_fpath) - - # fig = plt.figure(figsize=(9,6), dpi=108) - # ax1 = plt.subplot2grid((3,4), (0,0), colspan=2) - # ax1.set_title('EEG spectra before filtering') #TODO why it is not working? - # ax2 = plt.subplot2grid((3,4), (1,0), colspan=2) - # ax3 = plt.subplot2grid((3,4), (0,2), colspan=2) - # ax3.set_title('MEG spectra after filtering') - # ax4 = plt.subplot2grid((3,4), (1,2), colspan=2) - - # # raw.plot_psd(picks=['meg'], fmin=1, fmax=100, ax=[axes[0][0], axes[1][0]]) - # raw.plot_psd(picks=['meg'], fmin=1, fmax=100, ax=[ax1, ax2], show=False) - - # # raw_sss.plot_psd(picks=['meg'], fmin=1, fmax=100, ax=[axes[0][1], axes[1][1]]) - # raw_sss.plot_psd(picks=['meg'], fmin=1, fmax=100, ax=[ax3, ax4], show=False) - - # plt.axis('on') - # ax5 = plt.subplot2grid((3,4), (2,0), colspan=2) - # plt.axis('off') - # ax5.text(0, 0.7, 'noisy: ' + ', '.join(bad_chan['noisy'])) - # ax5.text(0, 0.4, 'flat: ' + ', '.join(bad_chan['flat'])) - - # pdf.savefig(fig) - # plt.close() - - ########################### - # Check EEG data quality # - ########################### - - if has_eeg: - print("has_eeg: viz_psd") - - fig = viz_psd(raw_sss) - pdf.savefig(fig) - plt.close() - - line_freqs = arange(raw.info['line_freq'], raw.info["sfreq"] / 2, raw.info['line_freq']) - - prep_params = { - "ref_chs": "eeg", - "reref_chs": "eeg", - "line_freqs": line_freqs, - "max_iterations": 4} - - montage = raw.get_montage() - prep = PrepPipeline(raw_sss, prep_params, montage, ransac=True) - prep.fit() - raw_car = prep.raw - raw_car.interpolate_bads(reset_bads=True) - - fig = viz_psd(raw_car) - pdf.savefig(fig) - plt.close() - print("end - has_eeg: viz_psd") - - #endregion PART 2a - MEG data filtering using Maxwell filters, method - defined in config - - #region annotations - - ########################### - # Detect ocular artifacts # - ########################### - - if has_eeg: - # Resetting the EOG channel - eog_ch = raw_sss.copy().pick_types(meg=False, eeg=False, eog=True) - if len(eog_ch.ch_names) < 2: - raw_sss.set_channel_types({'BIO002':'eog'}) - raw_sss.rename_channels({'BIO002': 'EOG002'}) - - # Find EOG events - eog_events = mne.preprocessing.find_eog_events(raw_sss) - onsets = (eog_events[:, 0] - raw_sss.first_samp) / raw_sss.info['sfreq'] - 0.25 - durations = [0.5] * len(eog_events) - descriptions = ['Blink'] * len(eog_events) - - # Annotate events - annot_blink = mne.Annotations( - onsets, - durations, - descriptions) - - ########################### - # Detect muscle artifacts # - ########################### - threshold_muscle = 7 - - # Notch filter - raw_muscle = raw_sss.copy().notch_filter([50, 100]) - - # Choose one channel type, if there are axial gradiometers and magnetometers, - # select magnetometers as they are more sensitive to muscle activity. - annot_muscle, scores_muscle = annotate_muscle_zscore( - raw_muscle, - ch_type="mag", - threshold=threshold_muscle, - min_length_good=0.3, - filter_freq=[110, 140]) - - ################# - # Detect breaks # - ################# - - # Get events - # events, event_id = mne.events_from_annotations(raw_sss) - - # Detect breaks based on events - # annot_break = mne.preprocessing.annotate_break( - # raw=raw_sss, - # events=events, - # min_break_duration=15.0) - - ########################### - - # Contatenate blink and muscle artifact annotations - if has_eeg: - annot_artifact = annot_blink + annot_muscle - else: - annot_artifact = annot_muscle - annot_artifact = mne.Annotations(onset = annot_artifact.onset + raw_sss._first_time, - duration = annot_artifact.duration, - description = annot_artifact.description, - orig_time = raw_sss.info['meas_date']) - - # Add artifact annotations in raw_sss - # raw_sss.set_annotations(raw_sss.annotations + annot_artifact + annot_break) - raw_sss.set_annotations(raw_sss.annotations + annot_artifact) - - #endregion annotations - - events, metadata = run_events(raw_sss, visit_id) - # Show events - fig = mne.viz.plot_events(events) - pdf.savefig(fig) - plt.close() - - raw_list.append(raw_sss) - events_list.append(events) - metadata_list.append(metadata) - - if run_part_3: - raw, events = mne.concatenate_raws(raw_list, events_list=events_list) - del raw_list - - # Concatenate metadata tables - metadata = pd.concat(metadata_list) - # metadata.to_csv(op.join(out_path, file_name[0:14] + 'ALL-meta.csv'), index=False) - - # Select sensor types - picks = mne.pick_types(raw.info, - meg = True, - eeg = has_eeg, - stim = True, - eog = has_eeg, - ecg = has_eeg, - ) - - # Set trial-onset event_ids - if visit_id == 'V1': - events_id = {} - types = ['face','object','letter','false'] - for j,t in enumerate(types): - for i in range(1,21): - events_id[t+str(i)] = i + j * 20 - elif visit_id == 'V2': - events_id = {} - events_id['blank'] = 50 - types = ['face','object'] - for j,t in enumerate(types): - for i in range(1,11): - events_id[t+str(i)] = i + j * 20 - - # Epoch raw data - epochs = mne.Epochs(raw, - events, - events_id, - tmin, tmax, - baseline=None, - proj=True, - picks=picks, - detrend=1, - reject=None, - reject_by_annotation=True, - verbose=True) - - del raw - - # Add metadata - epochs.metadata = metadata - - # Get rejection thresholds - reject = get_rejection_threshold(epochs, - ch_types=['mag', 'grad'], #'eeg'], #TODO: eeg not use for epoch rejection - decim=2) - - # Drop bad epochs based on peak-to-peak magnitude - epochs.drop_bad(reject=reject) - - # Plot percentage of rejected epochs per channel - fig1 = epochs.plot_drop_log() - pdf.savefig(fig1) - plt.close() - - -if __name__ == '__main__': - subject_id = input("Type the subject ID (e.g., SA101)\n>>> ") - visit_id = input("Type the visit ID (V1 or V2)\n>>> ") - run_qc_processing(subject_id, visit_id) \ No newline at end of file diff --git a/coglib/meeg/qc/srun_qc.py b/coglib/meeg/qc/srun_qc.py new file mode 100644 index 0000000..1dafd2b --- /dev/null +++ b/coglib/meeg/qc/srun_qc.py @@ -0,0 +1,29 @@ +#!/bin/bash +#SBATCH --partition=xnat +#SBATCH --nodes=1 +#SBATCH --cpus-per-task=2 +#SBATCH --mem-per-cpu=32G +#SBATCH --mail-type=BEGIN,END +#SBATCH --mail-user=gorska@wisc.edu +#SBATCH --time 12:00:00 +#SBATCH --chdir=/hpc/users/urszula.gorska/codes/MEEG/MNE-python_pipeline_v3/ + +if [ $# -ne 2 ]; + then echo "Please pass sub_prefix and visit as command line arguments. E.g." + echo "sbatch --array=101,103,105 srun_bids.sh SA V1" + echo "Exiting." + exit 1 +fi + +sub_prefix=$1 # Prefix of the subjects we're working on e.g. SA SB etc... +visit=$2 + +set -- + +module purge +module load Anaconda3/2020.11 +source /hpc/shared/EasyBuild/apps/Anaconda3/2020.11/bin/activate +conda activate /hpc/users/urszula.gorska/.conda/envs/mne_meg01_clone + +srun python P00_run_qc.py --sub ${sub_prefix}${SLURM_ARRAY_TASK_ID} --visit ${visit} +#srun python P00_run_qc_epochs.py --sub ${sub_prefix}`printf "%03d" $SLURM_ARRAY_TASK_ID` --visit ${visit} \ No newline at end of file From c7ca0f619c1218eb26cba39f26a49ae05a9410a0 Mon Sep 17 00:00:00 2001 From: Urszula Date: Wed, 5 Jul 2023 23:51:50 -0500 Subject: [PATCH 2/4] running qc --- coglib/meeg/README.md | 4 ++++ coglib/meeg/qc/srun_qc.py | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/coglib/meeg/README.md b/coglib/meeg/README.md index 993df63..031c41b 100644 --- a/coglib/meeg/README.md +++ b/coglib/meeg/README.md @@ -19,6 +19,10 @@ The environments are tailored for Linux and the HPC, so some things might break To run the analysis described below on the sample data, make sure to change the bids root path in /meeg/config/config.py: *$ROOT/sample_data/bids* +### Running 3rd level quality checks: +In the command line, enter: +sbatch --array= srun_qc.sh V1 + ### Running preprocessing: In the command line, enter: ``` diff --git a/coglib/meeg/qc/srun_qc.py b/coglib/meeg/qc/srun_qc.py index 1dafd2b..b64be28 100644 --- a/coglib/meeg/qc/srun_qc.py +++ b/coglib/meeg/qc/srun_qc.py @@ -16,7 +16,7 @@ fi sub_prefix=$1 # Prefix of the subjects we're working on e.g. SA SB etc... -visit=$2 +visit=$1 set -- From 3faf84c3aa8881d80f1daa301babb428ef33e571 Mon Sep 17 00:00:00 2001 From: Urszula Date: Thu, 6 Jul 2023 17:19:57 -0500 Subject: [PATCH 3/4] corrected file format --- coglib/meeg/qc/{srun_bids.py => srun_bids.sh} | 0 coglib/meeg/qc/{srun_qc.py => srun_qc.sh} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename coglib/meeg/qc/{srun_bids.py => srun_bids.sh} (100%) rename coglib/meeg/qc/{srun_qc.py => srun_qc.sh} (100%) diff --git a/coglib/meeg/qc/srun_bids.py b/coglib/meeg/qc/srun_bids.sh similarity index 100% rename from coglib/meeg/qc/srun_bids.py rename to coglib/meeg/qc/srun_bids.sh diff --git a/coglib/meeg/qc/srun_qc.py b/coglib/meeg/qc/srun_qc.sh similarity index 100% rename from coglib/meeg/qc/srun_qc.py rename to coglib/meeg/qc/srun_qc.sh From f2bb744b00df8c6d3a8634ddd845833c1b036bca Mon Sep 17 00:00:00 2001 From: Urszula Date: Thu, 6 Jul 2023 17:26:48 -0500 Subject: [PATCH 4/4] removing unnecessary files for sharing --- coglib/meeg/qc/QC_epochs.py | 198 ------------------------------------ coglib/meeg/qc/srun_bids.sh | 30 ------ coglib/meeg/qc/srun_qc.sh | 29 ------ 3 files changed, 257 deletions(-) delete mode 100644 coglib/meeg/qc/QC_epochs.py delete mode 100644 coglib/meeg/qc/srun_bids.sh delete mode 100644 coglib/meeg/qc/srun_qc.sh diff --git a/coglib/meeg/qc/QC_epochs.py b/coglib/meeg/qc/QC_epochs.py deleted file mode 100644 index 5ca3510..0000000 --- a/coglib/meeg/qc/QC_epochs.py +++ /dev/null @@ -1,198 +0,0 @@ -import os -import os.path as op -import matplotlib.pyplot as plt -import pandas as pd - - -import mne -import mne_bids -from pyprep.prep_pipeline import PrepPipeline -from mne.preprocessing import annotate_muscle_zscore -from autoreject import get_rejection_threshold -from numpy import arange - -from config import (bids_root, tmin, tmax) - -from matplotlib.backends.backend_pdf import PdfPages - -from qc.maxwell_filtering import run_maxwell_filter -from qc.extract_events import run_events - -from qc.viz_psd import viz_psd - - -def run_qc_epochs(subject_id, visit_id, has_eeg): - prep_deriv_root = op.join(bids_root, "derivatives", "preprocessing") - - # Set path to qc derivatives and create the related folders - qc_output_path = op.join(bids_root, "derivatives", "qc", visit_id) - if not op.exists(qc_output_path): - os.makedirs(qc_output_path) - - print("Processing subject: %s" % subject_id) - - #raw_list = list() - #events_list = list() - #metadata_list = list() - - # Set task - if visit_id == "V1": - bids_task = 'dur' - elif visit_id == "V2": - bids_task = 'vg' - elif visit_id == "V2": - bids_task = 'replay' - else: - raise ValueError("Error: could not set the task") - - - #with PdfPages(op.join(qc_output_path, subject_id + '_' + visit_id + '_MEG_V1_epochs.pdf')) as pdf: - - #FirstPage = plt.figure(figsize=(8,1), dpi=108) - #FirstPage.clf() - #plt.axis('off') - #plt.text(0.5, 0.5, subject_id, transform=FirstPage.transFigure, size=16, ha="center") - #pdf.savefig(FirstPage) - #plt.close() - - - #raw, events = mne.concatenate_raws(raw_list, events_list=events_list) - #del raw_list - - # Concatenate metadata tables - #metadata = pd.concat(metadata_list) - # metadata.to_csv(op.join(out_path, file_name[0:14] + 'ALL-meta.csv'), index=False) - - # Select sensor types - #picks = mne.pick_types(raw.info, - # meg = True, - # eeg = has_eeg, - # stim = True, - # eog = has_eeg, - # ecg = has_eeg, - # ) - - # Set trial-onset event_ids - if visit_id == 'V1': - events_id = {} - types = ['face','object','letter','false'] - for j,t in enumerate(types): - for i in range(1,21): - events_id[t+str(i)] = i + j * 20 -# elif visit_id == 'V2': -# events_id = {} -# events_id['blank'] = 50 -# types = ['face','object'] -# for j,t in enumerate(types): -# for i in range(1,11): -# events_id[t+str(i)] = i + j * 20 - - elif visit_id == 'V2': - if bids_task == "vg": - events_id = {} - events_id['blank'] = 50 - types = ['face','object'] - for j,t in enumerate(types): - for i in range(1,11): - events_id[t+str(i)] = i + j * 20 - elif bids_task == "replay": - events_id = {} - events_id['blankFT'] = 150 - events_id['blankOT'] = 250 - typesF = ['faceT','faceNT'] - for j,t in enumerate(typesF): - for i in range(1,11): - events_id[t+str(i)] = i + 100 + j * 100 - typesO = ['objectNT','objectT'] - for j,t in enumerate(typesO): - for i in range(1,11): - events_id[t+str(i)] = i + 120 + j * 100 - - # Epoch raw data - #epochs = mne.Epochs(raw, - # events, - # events_id, - # tmin, tmax, - # baseline=None, - # proj=True, - # picks=picks, - # detrend=1, - # reject=None, - # reject_by_annotation=True, - # verbose=True) - - # Add metadata - #epochs.metadata = metadata - - # Read epoched data from preprocessed - bids_path_epo = mne_bids.BIDSPath( - root=prep_deriv_root, - subject=subject_id, - datatype='meg', - task=bids_task, - session=visit_id, - suffix='epo', - extension='.fif', - check=False) - - epochs = mne.read_epochs(bids_path_epo.fpath, preload=False) - - if visit_id == 'V1': - print("VERY_IMPORTANT :)") - print("FACES task relevant") - epochs_rel_F = epochs['Task_relevance == "Relevant non-target" and Category == "face"'] - print(epochs_rel_F) - print("FACES task irrelevant") - epochs_irr_F = epochs['Task_relevance == "Irrelevant" and Category == "face"'] - print(epochs_irr_F) - - print("OBJECTS task relevant") - epochs_rel_O = epochs['Task_relevance == "Relevant non-target" and Category == "object"'] - print(epochs_rel_O) - print("OBJECTS task irrelevant") - epochs_irr_O = epochs['Task_relevance == "Irrelevant" and Category == "object"'] - print(epochs_irr_O) - - print("LETTERS task relevant") - epochs_rel_L = epochs['Task_relevance == "Relevant non-target" and Category == "letter"'] - print(epochs_rel_L) - print("LETTERS task irrelevant") - epochs_irr_L = epochs['Task_relevance == "Irrelevant" and Category == "letter"'] - print(epochs_irr_L) - - print("FALSE FONTS task relevant") - epochs_rel_S = epochs['Task_relevance == "Relevant non-target" and Category == "false"'] - print(epochs_rel_S) - print("FALSE FONTS task irrelevant") - epochs_irr_S = epochs['Task_relevance == "Irrelevant" and Category == "false"'] - print(epochs_irr_S) - - elif visit_id == 'V2': - print("FACES probe") - epochs_rel_F = epochs['Trial_type == "Probe" and Stimuli_type == "Face"'] - print(epochs_rel_F) - print("FACES filler") - epochs_irr_F = epochs['Trial_type == "Filler" and Stimuli_type == "Face"'] - print(epochs_irr_F) - - print("OBJECTS probe") - epochs_rel_O = epochs['Trial_type == "Probe" and Stimuli_type == "Object"'] - print(epochs_rel_O) - print("OBJECTS filler") - epochs_irr_O = epochs['Trial_type == "Filler" and Stimuli_type == "Object"'] - print(epochs_irr_O) - - print("BLANKS probe") - epochs_rel_O = epochs['Trial_type == "Probe" and Stimuli_type == "Blank"'] - print(epochs_rel_O) - print("BLANKS filler") - epochs_irr_O = epochs['Trial_type == "Filler" and Stimuli_type == "Blank"'] - print(epochs_irr_O) - - plt.close() - - -if __name__ == '__main__': - subject_id = input("Type the subject ID (e.g., SA101)\n>>> ") - visit_id = input("Type the visit ID (V1 or V2)\n>>> ") - run_qc_epochs(subject_id, visit_id) \ No newline at end of file diff --git a/coglib/meeg/qc/srun_bids.sh b/coglib/meeg/qc/srun_bids.sh deleted file mode 100644 index 688ac98..0000000 --- a/coglib/meeg/qc/srun_bids.sh +++ /dev/null @@ -1,30 +0,0 @@ -#!/bin/bash -#SBATCH --partition=xnat -#SBATCH --nodes=1 -#SBATCH --cpus-per-task=2 -#SBATCH --mem-per-cpu=32G -#SBATCH --mail-type=BEGIN,END -#SBATCH --mail-user=gorska@wisc.edu -#SBATCH --time 1:00:00 -#SBATCH --chdir=/hpc/users/urszula.gorska/codes/MEEG/MNE-python_pipeline_v3/ - -if [ $# -ne 2 ]; - then echo "Please pass sub_prefix visit and step as command line arguments. E.g." - echo "sbatch --array=101,103,105 srun_bids.sh SA V1" - echo "Exiting." - exit 1 -fi - -sub_prefix=$1 # Prefix of the subjects we're working on e.g. SA SB etc... -visit=$2 - -set -- - -module purge -module load Anaconda3/2020.11 -source /hpc/shared/EasyBuild/apps/Anaconda3/2020.11/bin/activate -conda activate /hpc/users/urszula.gorska/.conda/envs/mne_meg01_clone - -#srun python P00_bids_conversion.py --sub ${sub_prefix}${SLURM_ARRAY_TASK_ID} --visit ${visit} - -srun python P00_bids_conversion.py --sub ${sub_prefix}`printf "%03d" $SLURM_ARRAY_TASK_ID` --visit ${visit} \ No newline at end of file diff --git a/coglib/meeg/qc/srun_qc.sh b/coglib/meeg/qc/srun_qc.sh deleted file mode 100644 index b64be28..0000000 --- a/coglib/meeg/qc/srun_qc.sh +++ /dev/null @@ -1,29 +0,0 @@ -#!/bin/bash -#SBATCH --partition=xnat -#SBATCH --nodes=1 -#SBATCH --cpus-per-task=2 -#SBATCH --mem-per-cpu=32G -#SBATCH --mail-type=BEGIN,END -#SBATCH --mail-user=gorska@wisc.edu -#SBATCH --time 12:00:00 -#SBATCH --chdir=/hpc/users/urszula.gorska/codes/MEEG/MNE-python_pipeline_v3/ - -if [ $# -ne 2 ]; - then echo "Please pass sub_prefix and visit as command line arguments. E.g." - echo "sbatch --array=101,103,105 srun_bids.sh SA V1" - echo "Exiting." - exit 1 -fi - -sub_prefix=$1 # Prefix of the subjects we're working on e.g. SA SB etc... -visit=$1 - -set -- - -module purge -module load Anaconda3/2020.11 -source /hpc/shared/EasyBuild/apps/Anaconda3/2020.11/bin/activate -conda activate /hpc/users/urszula.gorska/.conda/envs/mne_meg01_clone - -srun python P00_run_qc.py --sub ${sub_prefix}${SLURM_ARRAY_TASK_ID} --visit ${visit} -#srun python P00_run_qc_epochs.py --sub ${sub_prefix}`printf "%03d" $SLURM_ARRAY_TASK_ID` --visit ${visit} \ No newline at end of file