diff --git a/analysis/1_physio.py b/analysis/1_physio.py index 8209cf4..2484544 100644 --- a/analysis/1_physio.py +++ b/analysis/1_physio.py @@ -34,7 +34,7 @@ def qc_physio(df, info, sub, plot_ecg=[], plot_ppg=[]): # Add text img = ill.image_text( - sub, color="black", size=100, x=-0.82, y=0.90, image=nk.fig2img(fig) + sub, color="black", size=40, x=-0.82, y=0.90, image=nk.fig2img(fig) ) plt.close(fig) # Do not show the plot in the console plot_ecg.append(img) @@ -51,7 +51,7 @@ def qc_physio(df, info, sub, plot_ecg=[], plot_ppg=[]): # Add text img = ill.image_text( - sub, color="black", size=100, x=-0.82, y=0.90, image=nk.fig2img(fig) + sub, color="black", size=40, x=-0.82, y=0.90, image=nk.fig2img(fig) ) plot_ppg.append(img) plt.close(fig) @@ -74,7 +74,7 @@ def qc_physio(df, info, sub, plot_ecg=[], plot_ppg=[]): # Loop through participants ================================================================== -for i, sub in enumerate(meta["participant_id"].values[84::]): +for i, sub in enumerate(meta["participant_id"].values): # Print progress and comments print(sub) @@ -93,104 +93,132 @@ def qc_physio(df, info, sub, plot_ecg=[], plot_ppg=[]): rs, hct = load_physio(path, sub) # Function loaded from script at URL # Resting State ========================================================================== + if sub not in ["sub-86"]: # No RS file + # Preprocessing -------------------------------------------------------------------------- + print(" - RS - Preprocessing") + + bio, info = nk.bio_process( + ecg=rs["ECG"][0][0], + ppg=rs["PPG_Muse"][0][0], + sampling_rate=rs.info["sfreq"], + ) + + # QC + qc["rs_ecg"], qc["rs_ppg"] = qc_physio( + bio, info, sub, plot_ecg=qc["rs_ecg"], plot_ppg=qc["rs_ppg"] + ) + + # Hear Rate Variability (HRV) ------------------------------------------------------------- + hrv = nk.hrv( + bio["ECG_R_Peaks"].values.nonzero()[0], sampling_rate=rs.info["sfreq"] + ) + idx = [ + "MeanNN", + "SDNN", + "RMSSD", + "SampEn", + "HF", + "HFD", + "LFHF", + "IALS", + "Ca", + "AI", + ] + hrv = hrv[["HRV_" + s for s in idx]] + dfsub = pd.concat([dfsub, hrv], axis=1) + + # Heartbeat Counting Task (HCT) =========================================================== # Preprocessing -------------------------------------------------------------------------- - print(" - RS - Preprocessing") + print(" - HCT - Preprocessing") + + # Find events (again as data was cropped) and epoch + events = nk.events_find( + hct["PHOTO"][0][0], threshold_keep="below", duration_min=15000 + ) - bio, info = nk.bio_process( - ecg=rs["ECG"][0][0], - ppg=rs["PPG_Muse"][0][0], - sampling_rate=rs.info["sfreq"], + # Make sure there are 6 events of expected duration + # - sub68 has one slitghtly shorter event (we can keep it) + if sub not in ["sub-68"]: + assert len(events["onset"]) == 6 + durations = np.sort(events["duration"] / hct.info["sfreq"]) + assert np.max(np.abs(durations - np.array([20, 25, 30, 35, 40, 45]))) < 0.75 + + # Process signals + hct_physio = hct.to_data_frame()[["PHOTO", "ECG", "PPG_Muse"]] + hct_physio, info = nk.bio_process( + ecg=hct_physio["ECG"].values, + ppg=hct_physio["PPG_Muse"].values, + sampling_rate=hct.info["sfreq"], + keep=hct_physio["PHOTO"], ) # QC - qc["rs_ecg"], qc["rs_ppg"] = qc_physio( - bio, info, sub, plot_ecg=qc["rs_ecg"], plot_ppg=qc["rs_ppg"] + qc["hct_ecg"], qc["hct_ppg"] = qc_physio( + hct_physio, info, sub, plot_ecg=qc["hct_ecg"], plot_ppg=qc["hct_ppg"] ) + # Analysis -------------------------------------------------------------------------- + # Make epochs + epochs = nk.epochs_create( + hct_physio, + events, + sampling_rate=hct.info["sfreq"], + epochs_start=0, + epochs_end="from_events", + ) + + # Load behavioral data + file = [file for file in os.listdir(path_beh) if "HCT" in file] + file = path_beh + [f for f in file if ".tsv" in f][0] + hct_beh = pd.read_csv(file, sep="\t") + + # Count R peaks in each epoch + hct_beh["N_R_peaks"] = [epoch["ECG_R_Peaks"].sum() for i, epoch in epochs.items()] + hct_beh["N_PPG_peaks"] = [epoch["PPG_Peaks"].sum() for i, epoch in epochs.items()] + + # Manual fix (based on comments) + peaks = hct_beh["N_R_peaks"].values + if sub == "sub-09": # No ECG + peaks = hct_beh["N_PPG_peaks"].values + + # Compute accuracy + if sub in ["sub-72"]: + hct_beh["HCT_Accuracy"] = np.nan + else: + hct_beh["HCT_Accuracy"] = 1 - ((np.abs(hct_beh["Answer"] - peaks)) / peaks) + + # Manual fixes (based on comments) + if sub == "sub-07": + hct_beh.loc[2, ["Confidence", "HCT_Accuracy"]] = np.nan + if sub == "sub-11": + hct_beh.loc[0:2, ["Confidence", "HCT_Accuracy"]] = np.nan + + # Compute interoception scores (Garfinkel et al., 2015) ----------------------------------- + if sub not in ["sub-72"]: + dfsub["HCT_Accuracy"] = np.nanmean(hct_beh["HCT_Accuracy"]) + dfsub["HCT_Sensibility"] = np.nanmean(hct_beh["Confidence"]) + dfsub["HCT_Awareness"] = scipy.stats.spearmanr( + hct_beh["Confidence"].dropna(), hct_beh["HCT_Accuracy"].dropna() + ).statistic + # Hear Rate Variability (HRV) ------------------------------------------------------------- - hrv = nk.hrv(bio["ECG_R_Peaks"].values.nonzero()[0], sampling_rate=rs.info["sfreq"]) - idx = ["MeanNN", "SDNN", "RMSSD", "SampEn", "HF", "HFD", "LFHF", "IALS", "Ca", "AI"] - hrv = hrv[["HRV_" + s for s in idx]] - dfsub = pd.concat([dfsub, hrv], axis=1) - - # # Heartbeat Counting Task (HCT) =========================================================== - - # # Preprocessing -------------------------------------------------------------------------- - # print(" - HCT - Preprocessing") - - # # Find events (again as data was cropped) and epoch - # events = nk.events_find( - # hct["PHOTO"][0][0], threshold_keep="below", duration_min=15000 - # ) - - # # Process signals - # hct_physio = hct.to_data_frame()[["PHOTO", "ECG", "PPG_Muse"]] - # hct_physio, info = nk.bio_process( - # ecg=hct_physio["ECG"].values, - # ppg=hct_physio["PPG_Muse"].values, - # sampling_rate=2000, - # keep=hct_physio["PHOTO"], - # ) - - # # # QC - # # qc["hct_ecg"], qc["hct_ppg"] = qc_physio( - # # hct_physio, info, sub, plot_ecg=qc["hct_ecg"], plot_ppg=qc["hct_ppg"] - # # ) - - # # Analysis -------------------------------------------------------------------------- - # # Make epochs - # epochs = nk.epochs_create( - # hct_physio, events, sampling_rate=2000, epochs_start=0, epochs_end="from_events" - # ) - - # # Load behavioral data - # file = [file for file in os.listdir(path_beh) if "HCT" in file] - # file = path_beh + [f for f in file if ".tsv" in f][0] - # hct_beh = pd.read_csv(file, sep="\t") - - # # Count R peaks in each epoch - # hct_beh["N_R_peaks"] = [epoch["ECG_R_Peaks"].sum() for i, epoch in epochs.items()] - # hct_beh["N_PPG_peaks"] = [epoch["PPG_Peaks"].sum() for i, epoch in epochs.items()] - - # # Manual fix (based on comments) - # peaks = hct_beh["N_R_peaks"].values - # if sub == "sub-09": # No ECG - # peaks = hct_beh["N_PPG_peaks"].values - - # # Compute accuracy - # hct_beh["HCT_Accuracy"] = 1 - ((np.abs(hct_beh["Answer"] - peaks)) / peaks) - - # # Manual fixes (based on comments) - # if sub == "sub-07": - # hct_beh.loc[2, ["Confidence", "HCT_Accuracy"]] = np.nan - # if sub == "sub-11": - # hct_beh.loc[0:2, ["Confidence", "HCT_Accuracy"]] = np.nan - - # # Compute interoception scores (Garfinkel et al., 2015) ----------------------------------- - # dfsub["HCT_Accuracy"] = np.mean(hct_beh["HCT_Accuracy"]) - # dfsub["HCT_Sensibility"] = np.nanmean(hct_beh["Confidence"]) - # dfsub["HCT_Awareness"] = scipy.stats.spearmanr( - # hct_beh["Confidence"].dropna(), hct_beh["HCT_Accuracy"].dropna() - # ).statistic - - # # Hear Rate Variability (HRV) ------------------------------------------------------------- - # print(" - HCT - HRV") - - # hrv = pd.concat( - # [ - # nk.hrv_time( - # e["ECG_R_Peaks"].values.nonzero()[0], sampling_rate=rs.info["sfreq"] - # ) - # for e in epochs.values() - # ] - # ) - # hrv = hrv[["HRV_" + s for s in ["MeanNN", "SDNN", "RMSSD"]]].mean(axis=0) - # hrv.index = [s + "_HCT" for s in hrv.index] - # dfsub = pd.concat([dfsub, pd.DataFrame(hrv).T], axis=1) - - # # Append participant to rest -------------------------------------------------------------- - # df = pd.concat([df, dfsub], axis=0) + print(" - HCT - HRV") + + hrv = pd.concat( + [ + nk.hrv_time( + e["ECG_R_Peaks"].values.nonzero()[0], sampling_rate=hct.info["sfreq"] + ) + for e in epochs.values() + ] + ) + hrv = hrv[["HRV_" + s for s in ["MeanNN", "SDNN", "RMSSD"]]].mean(axis=0) + hrv.index = [s + "_HCT" for s in hrv.index] + dfsub = pd.concat([dfsub, pd.DataFrame(hrv).T], axis=1) + + # Append participant to rest -------------------------------------------------------------- + df = pd.concat([df, dfsub], axis=0) # Clean up and Save data diff --git a/analysis/func_load_physio.py b/analysis/func_load_physio.py index 83ccb00..6e4678a 100644 --- a/analysis/func_load_physio.py +++ b/analysis/func_load_physio.py @@ -21,73 +21,83 @@ def consecutive_nans(raw): other[ch] = consecutive return np.max(start), other - # Path to EEG data - path_eeg = path + sub + "/eeg/" - path_beh = path + sub + "/beh/" - - # Resting State ============================================================================== - file = [file for file in os.listdir(path_eeg) if "RS" in file] - file = path_eeg + [f for f in file if ".vhdr" in f][0] - - rs = mne.io.read_raw_brainvision(file, preload=True) - rs = rs.set_montage("standard_1020") - - # Detect onset of RS - events = nk.events_find( - rs.to_data_frame()["PHOTO"], # nk.signal_plot(rs["PHOTO"][0][0]) - threshold_keep="below", - duration_min=int(rs.info["sfreq"] * 5), - ) + # RS ============================================================================== + def load_rs(path, sub): + + # Path to EEG data + path_eeg = path + sub + "/eeg/" + file = [file for file in os.listdir(path_eeg) if "RS" in file] + file = path_eeg + [f for f in file if ".vhdr" in f][0] + + rs = mne.io.read_raw_brainvision(file, preload=True) + rs = rs.set_montage("standard_1020") + + # Detect onset of RS + events = nk.events_find( + rs.to_data_frame()["PHOTO"], # nk.signal_plot(rs["PHOTO"][0][0]) + threshold_keep="below", + duration_min=int(rs.info["sfreq"] * 5), + ) - if sub in ["sub-57"]: # No Lux, crop out based on GYRO - # nk.signal_plot(rs["GYRO"][0][0]) - # plt.vlines( - # [10000, 10000 + rs.info["sfreq"] * 60 * 8], ymin=0, ymax=10, color="red" - # ) - events = {"onset": [10000], "duration": [rs.info["sfreq"] * 60 * 8]} - if sub in ["sub-80"]: # No Lux, crop out based on GYRO + # No Lux, crop out based on GYRO # nk.signal_plot(rs["GYRO"][0][0]) # plt.vlines( # [10000, 10000 + rs.info["sfreq"] * 60 * 8], ymin=0, ymax=10, color="red" # ) - events = {"onset": [10000], "duration": [rs.info["sfreq"] * 60 * 8]} - assert len(events["onset"]) == 1 # Check that there is only one event - - rs = nk.mne_crop( - rs, smin=events["onset"][0], smax=events["onset"][0] + events["duration"][0] - ) - assert len(rs) / rs.info["sfreq"] / 60 > 6 # Check duration is at least 6 minutes - - # Cut out nans at the beginning due to sync delay - if sub in ["sub-24", "sub-38", "sub-65", "sub-68", "sub-76", "sub-105"]: - # rs.to_data_frame() - first_valid, _ = consecutive_nans(rs) - rs = nk.mne_crop(rs, smin=first_valid, smax=None) - assert rs.to_data_frame()["ECG"].isna().sum() == 0 - - # Check MUSE signal interruptions - _, others = consecutive_nans(rs) - if sub in [ - "sub-10", - "sub-15", - "sub-16", - "sub-19", - "sub-20", - "sub-31", - "sub-42", - "sub-50", - "sub-70", - "sub-76", - ]: - rs = rs.apply_function( - nk.signal_fillmissing, picks="PPG_Muse", method="forward" + if sub in ["sub-57", "sub-80", "sub-91"]: + events = {"onset": [10000], "duration": [rs.info["sfreq"] * 60 * 8]} + if sub in ["sub-93"]: + events = {"onset": [4600], "duration": [rs.info["sfreq"] * 60 * 8]} + assert len(events["onset"]) == 1 # Check that there is only one event + + rs = nk.mne_crop( + rs, smin=events["onset"][0], smax=events["onset"][0] + events["duration"][0] ) + assert ( + len(rs) / rs.info["sfreq"] / 60 > 6 + ) # Check duration is at least 6 minutes + + # Cut out nans at the beginning due to sync delay + if sub in ["sub-24", "sub-38", "sub-65", "sub-68", "sub-76", "sub-105"]: + # rs.to_data_frame() + first_valid, _ = consecutive_nans(rs) + rs = nk.mne_crop(rs, smin=first_valid, smax=None) + assert rs.to_data_frame()["ECG"].isna().sum() == 0 + + # Check MUSE signal interruptions + _, others = consecutive_nans(rs) + if sub in [ + "sub-10", + "sub-15", + "sub-16", + "sub-19", + "sub-20", + "sub-31", + "sub-42", + "sub-50", + "sub-70", + "sub-76", + "sub-82", + "sub-95", + ]: + rs = rs.apply_function( + nk.signal_fillmissing, picks="PPG_Muse", method="forward" + ) + else: + # rs.to_data_frame(["ECG", "AF7", "PPG_Muse", "PHOTO"]).plot(subplots=True) + assert ( + len(others["PPG_Muse"]) == 0 + ) # Make sure no missing values for PPG-Muse + return rs + + if sub in ["sub-86"]: # No RS file + rs = None else: - # rs.to_data_frame(["ECG", "AF7", "PPG_Muse", "PHOTO"]).plot(subplots=True) - assert len(others["PPG_Muse"]) == 0 # Make sure no missing values for PPG-Muse + rs = load_rs(path, sub) # HCT ============================================================================== # Open HCT file + path_eeg = path + sub + "/eeg/" file = [file for file in os.listdir(path_eeg) if "HCT" in file] file = path_eeg + [f for f in file if ".vhdr" in f][0] hct = mne.io.read_raw_brainvision(file, preload=True, verbose=False) @@ -125,7 +135,17 @@ def consecutive_nans(raw): # Check if MUSE signal interruptions _, others = consecutive_nans(hct) - if sub in ["sub-03", "sub-04", "sub-11", "sub-12", "sub-13", "sub-70", "sub-103"]: + if sub in [ + "sub-03", + "sub-04", + "sub-11", + "sub-12", + "sub-13", + "sub-70", + "sub-89", + "sub-95", + "sub-103", + ]: hct = hct.apply_function( nk.signal_fillmissing, picks="PPG_Muse", method="forward" )