Merge pull request #119 from kymata-atlas/big-preprocessing-refactor-cov

Add the function to calculate noise covariance matrix

invokers/invoker_estimate_noise_cov_matrix.py

@@ -0,0 +1,31 @@
+from pathlib import Path
+from colorama import Fore
+
+from kymata.io.yaml import load_config
+from kymata.io.cli import print_with_color
+from kymata.preproc.data_cleansing import estimate_noise_cov
+
+
+# noinspection DuplicatedCode
+def main():
+    config = load_config(str(Path(Path(__file__).parent.parent, "kymata", "config", "dataset4.yaml")))
+
+    if config['data_location'] == "local":
+        data_root_dir = str(Path(Path(__file__).parent.parent, "kymata-toolbox-data", "emeg_study_data")) + "/"
+    elif config['data_location'] == "cbu":
+        data_root_dir = '/imaging/projects/cbu/kymata/data/'
+    elif config['data_location'] == "cbu-local":
+        data_root_dir = '//cbsu/data/imaging/projects/cbu/kymata/data/'
+    else:
+        raise Exception("The 'data_location' parameter in the config file must be either 'cbu' or 'local' or 'cbu-local'.")
+
+    estimate_noise_cov( data_root_dir = data_root_dir,
+                       emeg_machine_used_to_record_data = config['emeg_machine_used_to_record_data'],
+                       list_of_participants = config['list_of_participants'],
+                       dataset_directory_name  = config['dataset_directory_name'],
+                       n_runs = config['number_of_runs'],
+                       cov_method = config['cov_method']
+                        )
+
+if __name__ == '__main__':
+    main()

invokers/invoker_run_data_cleansing.py

@@ -17,7 +17,7 @@ def main():
     elif config['data_location'] == "cbu-local":
         data_root_dir = '//cbsu/data/imaging/projects/cbu/kymata/data/'
     else:
-        raise Exception("The 'data_location' parameter in the config file must be either 'cbu' or 'local'.")
+        raise Exception("The 'data_location' parameter in the config file must be either 'cbu' or 'local' or 'cbu-local'.")
 
     run_first_pass_cleansing_and_maxwell_filtering(
         data_root_dir = data_root_dir,

invokers/invoker_run_hexel_current_estimation.py

@@ -17,9 +17,9 @@ def main():
     elif config['data_location'] == "cbu-local":
         data_root_dir = '//cbsu/data/imaging/projects/cbu/kymata/data/'
     else:
-        raise Exception("The 'data_location' parameter in the config file must be either 'cbu' or 'local'.")
+        raise Exception("The 'data_location' parameter in the config file must be either 'cbu' or 'local' or 'cbu-local'.")
 
-    create_current_estimation_prerequisites(data_root_dir, config=config)
+    # create_current_estimation_prerequisites(data_root_dir, config=config)
 
     create_forward_model_and_inverse_solution(data_root_dir, config=config)
 

kymata/config/dataset4.yaml

@@ -32,6 +32,9 @@ supress_excessive_plots_and_prompts: True
 eeg: True
 meg: True
 
+# Method to estimate noise covariance matrix
+cov_method: 'grand_ave'            # grand_ave | empty_room | run_start
+
 # Creation of evoked trial pipeline
 audio_delivery_latency:  26        # in milliseconds
 visual_delivery_latency: 17        # in milliseconds

kymata/preproc/data_cleansing.py

@@ -5,6 +5,7 @@
 from numpy import zeros, nanmin
 from pandas import DataFrame, Index
 from pathlib import Path
+from matplotlib import pyplot as plt
 
 from kymata.io.cli import print_with_color, input_with_color
 from kymata.io.yaml import load_config, modify_param_config
@@ -297,6 +298,43 @@ def _remove_ecg_eog(filt_raw, ica):
     ica.plot_scores(eog_scores)
 
 
+def estimate_noise_cov(data_root_dir: str,
+                       emeg_machine_used_to_record_data: str,
+                       list_of_participants: list[str],
+                       dataset_directory_name: str,
+                       n_runs: int,
+                       cov_method: str,
+                       ):
+    for p in list_of_participants:
+        if cov_method == 'grand_ave':
+            cleaned_raws = []
+            for run in range(1, n_runs + 1):
+                raw_fname = data_root_dir + dataset_directory_name + '/intrim_preprocessing_files/2_cleaned/' + p + '_run' + str(run) + '_cleaned_raw.fif.gz'
+                raw = mne.io.Raw(raw_fname, preload=True)
+                cleaned_raws.append(raw)
+            raw = mne.io.concatenate_raws(raws=cleaned_raws, preload=True)
+            cov = mne.compute_raw_covariance(raw, tmin=300, tmax=310, return_estimators=True)
+            mne.write_cov(data_root_dir + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/covariance_grand_average/' + p + '-auto-cov-300.fif', cov)
+
+        elif cov_method == 'empty_room':
+            emptyroom_raw = mne.read_raw() 
+            emptyroom_raw = mne.preprocessing.maxwell_filter_prepare_emptyroom(emptyroom_raw)
+
+            fine_cal_file = str(Path(Path(__file__).parent.parent.parent, 'kymata-toolbox-data', 'cbu_specific_files/SSS/sss_cal_' + emeg_machine_used_to_record_data + '.dat'))
+            crosstalk_file = str(Path(Path(__file__).parent.parent.parent, 'kymata-toolbox-data', 'cbu_specific_files/SSS/ct_sparse_' + emeg_machine_used_to_record_data + '.fif'))
+
+
+            raw_fif_data_sss = mne.preprocessing.maxwell_filter(
+                        emptyroom_raw,
+                        cross_talk=crosstalk_file,
+                        calibration=fine_cal_file,
+                        st_correlation=0.980,
+                        st_duration=10,
+                        verbose=True)
+
+            cov = mne.compute_raw_covariance(raw_fif_data_sss, tmin=0, tmax=10, return_estimators=True)
+            mne.write_cov(data_root_dir + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/covariance_grand_average/' + p + '-auto-cov-emptyroom.fif', cov)
+
 def create_trialwise_data(dataset_directory_name: str,
                         list_of_participants: list[str],
                         repetitions_per_runs: int,
@@ -414,8 +452,8 @@ def create_trialwise_data(dataset_directory_name: str,
             evoked = epochs.average()
             evoked.save(f'{data_path}/intrim_preprocessing_files/{save_folder}/evoked_data/{p}-ave.fif', overwrite=True)
 
-
-def create_trials(dataset_directory_name: str,
+def create_trials(data_root_dir: str,
+                  dataset_directory_name: str,
                   list_of_participants: list[str],
                   repetitions_per_runs: int,
                   number_of_runs: int,
@@ -443,7 +481,7 @@ def create_trials(dataset_directory_name: str,
         cleaned_raws = []
 
         for run in range(1, number_of_runs + 1):
-            raw_fname = '/imaging/projects/cbu/kymata/data/' + dataset_directory_name + '/intrim_preprocessing_files/2_cleaned/' + p + '_run' + str(run) + '_cleaned_raw.fif.gz'
+            raw_fname = data_root_dir + dataset_directory_name + '/intrim_preprocessing_files/2_cleaned/' + p + '_run' + str(run) + '_cleaned_raw.fif.gz'
             raw = mne.io.Raw(raw_fname, preload=True)
             cleaned_raws.append(raw)
 
@@ -510,7 +548,7 @@ def create_trials(dataset_directory_name: str,
             # 	Log which channels are worst
             dropfig = epochs.plot_drop_log(subject=p)
             dropfig.savefig(
-                '/imaging/projects/cbu/kymata/data/' + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/logs/' + input_stream + '_drop-log_' + p + '.jpg')
+                data_root_dir + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/logs/' + input_stream + '_drop-log_' + p + '.jpg')
 
             global_droplog.append('[' + input_stream + ']' + p + ':' + str(epochs.drop_log_stats(epochs.drop_log)))
 
@@ -532,13 +570,6 @@ def create_trials(dataset_directory_name: str,
 #            'data/' + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/evoked_grand_average/' + p + '-grandave.fif',
 #            overwrite=True)
 
-        # save grand covs
-        print_with_color(f"... save grand covariance matrix", Fore.GREEN)
-
-        cov = mne.compute_raw_covariance(raw, tmin=0, tmax=10, return_estimators=True)
-        mne.write_cov('/imaging/projects/cbu/kymata/data/' + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/covariance_grand_average/' + p + '-auto-cov.fif', cov)
-
-
 # Save global droplog
 #    with open('data/' + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/logs/drop-log.txt', 'a') as file:
 #        file.write('Average drop rate for each participant\n')
@@ -563,6 +594,7 @@ def create_trials(dataset_directory_name: str,
 #                'data/' + dataset_directory_name + '/intrim_preprocessing_files/3_evoked_sensor_data/evoked_data/' + input_stream + '/item' + str(
 #                    trial) + '-ave.fif', overwrite=True)
 
+
 def apply_automatic_bad_channel_detection(raw_fif_data: mne.io.Raw, machine_used: str, plot: bool = True):
     """Apply Automatic Bad Channel Detection"""
     raw_check = raw_fif_data.copy()
@@ -588,10 +620,7 @@ def apply_automatic_bad_channel_detection(raw_fif_data: mne.io.Raw, machine_used
 
     return raw_fif_data
 
-
 def _plot_bad_chans(auto_scores):
-    import seaborn as sns
-    from matplotlib import pyplot as plt
 
     # Only select the data for gradiometer channels.
     ch_type = 'grad'
@@ -635,4 +664,4 @@ def _plot_bad_chans(auto_scores):
     fig.tight_layout(rect=[0, 0.03, 1, 0.95])
 
     # Replace the word “noisy” with “flat”, and replace
-    # vmin=nanmin(limits) with vmax=nanmax(limits) to print flat channels
+    # vmin=nanmin(limits) with vmax=nanmax(limits) to print flat channels

kymata/preproc/hexel_current_estimation.py

@@ -171,33 +171,33 @@ def create_forward_model_and_inverse_solution(data_root_dir, config: dict):
     intrim_preprocessing_directory_name = Path(data_root_dir, dataset_directory_name, "intrim_preprocessing_files")
     mri_structurals_directory = Path(data_root_dir, dataset_directory_name, config['mri_structurals_directory'])
 
-    # Compute forward solution
-    for participant in list_of_participants:
+    # # Compute forward solution
+    # for participant in list_of_participants:
 
-        fwd = mne.make_forward_solution(
-            # Path(Path(path.abspath("")), "data",
-            Path(data_root_dir,
-                 dataset_directory_name,
-                 'raw_emeg', participant, participant +
-                 '_run1_raw.fif'), # note this file is only used for the sensor positions.
-            trans=Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","coregistration_files", participant + '-trans.fif'),
-            src=Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","src_files", participant + '_ico5-src.fif'),
-            bem=Path(mri_structurals_directory, participant, "bem", participant + '-5120-5120-5120-bem-sol.fif'),
-            meg=config['meg'],
-            eeg=config['eeg'],
-            mindist=5.0,
-            n_jobs=None,
-            verbose=True,
-        )
-        print(fwd)
-        if config['meg'] and config['eeg']:
-            mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd.fif'), fwd)
-        elif config['meg']:
-            mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd-megonly.fif'), fwd)
-        elif config['eeg']:
-            mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd-eegonly.fif'), fwd)
-        else:
-            raise Exception('eeg and meg in the config file cannot be both False')
+    #     fwd = mne.make_forward_solution(
+    #         # Path(Path(path.abspath("")), "data",
+    #         Path(data_root_dir,
+    #              dataset_directory_name,
+    #              'raw_emeg', participant, participant +
+    #              '_run1_raw.fif'), # note this file is only used for the sensor positions.
+    #         trans=Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","coregistration_files", participant + '-trans.fif'),
+    #         src=Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","src_files", participant + '_ico5-src.fif'),
+    #         bem=Path(mri_structurals_directory, participant, "bem", participant + '-5120-5120-5120-bem-sol.fif'),
+    #         meg=config['meg'],
+    #         eeg=config['eeg'],
+    #         mindist=5.0,
+    #         n_jobs=None,
+    #         verbose=True,
+    #     )
+    #     print(fwd)
+    #     if config['meg'] and config['eeg']:
+    #         mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd.fif'), fwd)
+    #     elif config['meg']:
+    #         mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd-megonly.fif'), fwd)
+    #     elif config['eeg']:
+    #         mne.write_forward_solution(Path(intrim_preprocessing_directory_name, "4_hexel_current_reconstruction","forward_sol_files", participant + '-fwd-eegonly.fif'), fwd)
+    #     else:
+    #         raise Exception('eeg and meg in the config file cannot be both False')
 
     # Compute inverse operator
 
@@ -224,7 +224,7 @@ def create_forward_model_and_inverse_solution(data_root_dir, config: dict):
             intrim_preprocessing_directory_name,
             '3_evoked_sensor_data',
             'covariance_grand_average',
-            participant + '-auto-cov.fif')))
+            participant + '-auto-cov-300.fif')))
         # note this file is only used for the sensor positions.
         raw = mne.io.Raw(Path(
             Path(path.abspath("")),
@@ -249,7 +249,7 @@ def create_forward_model_and_inverse_solution(data_root_dir, config: dict):
                     intrim_preprocessing_directory_name,
                     '4_hexel_current_reconstruction',
                     'inverse-operators',
-                    participant + '_ico5-3L-loose02-cps-nodepth.fif')), 
+                    participant + '_ico5-3L-loose02-cps-nodepth-300.fif')), 
                 inverse_operator)
         elif config['meg']:
             mne.minimum_norm.write_inverse_operator(