<jupyter_start><jupyter_text>Tutorial 3: Regressors from ecg data This example shows how to use *niphlem* for one subject from one of our particular projects. These data consist of one acquisition of resting-state that has been previously preprocesses using *fMRIPrep* and physiological acquired through a Siemens multiband protocol. For the sake of simplicity in this example, we are going to concentrate only on cardiac signal reprsented by the ECG. In addition, we are going to show the maps generated by the *niphlem* generated regressors using the library *nilearn* (https://nilearn.github.io/)<jupyter_code># IMPORTS
import numpy as np
import pandas as pd

from nilearn import image, plotting
from nilearn.glm.first_level import FirstLevelModel

from niphlem.input_data import (load_cmrr_info, load_cmrr_data)<jupyter_output><empty_output><jupyter_text>Let's start definining our data<jupyter_code># Our preprocessed BOLD image
run_img = image.load_img("./data/demo/sub-06_ses-04_task-resting_run-01_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz")
# The repetition time of the scanner
t_r = 1.5
# The frame times, used to define the regressors at the precise times
n_vols = run_img.shape[-1]
frame_times = np.arange(n_vols)*t_r
# Our file with the confounder regressors as a Dataframe. It will be used later to extract the motion parameters
conf_df = pd.read_csv("./data/demo/sub-06_ses-04_task-resting_run-01_desc-confounds_regressors.tsv", sep="\t")

# Our physiological data, the info log file and the ECG log file
info_log = "./data/demo/physio/Physio_20210322_140315_89a222d1-4c24-4caf-a898-f06c6bfd2342_Info.log"
ecg_log = "./data/demo/physio/Physio_20210322_140315_89a222d1-4c24-4caf-a898-f06c6bfd2342_ECG.log"<jupyter_output><empty_output><jupyter_text>Load data The first thing that we need is to load the info file, which will give us the times of the scanner. We can load this kind of the data using the function `niphlem.input_data.load_cmrr_info`. The only thing that we have to do is to pass the info log file to this function, and it will return the time traces and a dictionary with some meta information which is at the beginning of the file.<jupyter_code>time_traces, meta_info = load_cmrr_info(info_log)
meta_info<jupyter_output><empty_output><jupyter_text>Now that we have the time traces, we are going to load our ECG data by passing to the function `niphlem.input_data.load_cmrr_data` its log file and the meta information dictionary that we've just generated. This dictionary is required to be able to synchronize the physiological signal with the scanner times. We also specify the type of signal, as a new entry to an updated dictionary will be added, containing meta information of this ECG signal.<jupyter_code>ecg_signal, meta_info = load_cmrr_data(ecg_log, info_dict=meta_info, sig_type="cardiac")<jupyter_output><empty_output><jupyter_text>Let's have a look into a portion of this signal, corresponding to 2 s<jupyter_code>fig, ax =plt.subplots(figsize=(7,5))
plt.plot(ecg_signal[:,0], label="electrode 1")
plt.plot(1000 + ecg_signal[:,1], label="electrode 2")
plt.plot(2000 + ecg_signal[:,2], label="electrode 3")
plt.plot(3000 + ecg_signal[:,3], label="electrode 4")
plt.legend(prop={'size':12})
plt.xlim([0,800])
pass<jupyter_output><empty_output><jupyter_text>Let's also have a look at the updated information dictionary.<jupyter_code>meta_info<jupyter_output><empty_output><jupyter_text>Retroicor regressors Let's start the generation of regressors for the Retroicor Model. The class for generating these regressors is found as `niphlem.models.RetroicorPhysio`.<jupyter_code>from niphlem.models import RetroicorPhysio<jupyter_output><empty_output><jupyter_text>Let's have a look at the documentation of the class<jupyter_code>print(RetroicorPhysio.__doc__)<jupyter_output>Physiological regressors using Retroicor.

    Parameters
    ----------
    physio_rate : float
        Sampling rate for the recording in Hz.
        This is needed for filtering to define the nyquist frequency.
    t_r : float
        Repetition time for the scanner (the usual T_R) in secs.
    delta: float
        minimum separation (in physio recording units) between
        events in signal to be considered peaks
    peak_rise: float
        relative height with respect to the 20th tallest events in signal
        to consider events as peak.
    order: int or array-like of shape (n_orders,)
        Fourier expansion for phases. If int, the fourier expansion is
        performed to that order, starting from 1. If an array is provided,
        each element will multiply the phases.
    transform : {"demean", "zscore", "abs"}, optional
        Transform data before filtering. The default is "demean".
    high_pass : float, optional
        High-pass filtering frequency (in Hz)[...]<jupyter_text>We have three required arguments (physio_rate, t_r, delta) and a bunch of optional ones. In our case *physio_rate* is 400 Hz, t_r is 1.5 s and delta=200, the minimu separation between peaks and that corresponds to 200, i.e. 2 peaks per second or at most 120 peaks per minute. Before the computing any regressors, we are going to take the mean signal across all the for electrodes and then filter the signal between 0.5 and 3.0 Hz.Finally, retroicor regressors of order 2 will be generated.<jupyter_code>retro_ecg = RetroicorPhysio(physio_rate=400, # Physiological rate (400 Hz)
                            t_r=1.5, # Scan rate (1.5 s)
                            delta=200, 
                            peak_rise=0.5,
                            columns="mean", # mean across the four electrodes
                            low_pass=2.0, # keep frequencies below this
                            high_pass=0.5, # keep frequencies above this
                            order=2 # order 2 of the retroicor expansion
                           )<jupyter_output><empty_output><jupyter_text>Once we have defined this object, we generate the regressors using the method `compute_regressors`, which is comon to all the models<jupyter_code>retroicor_regressors = retro_ecg.compute_regressors(signal=ecg_signal, time_scan=frame_times)<jupyter_output><empty_output><jupyter_text>Let's see the a part of these regressors<jupyter_code>fig, ax = plt.subplots(figsize=(7,5))
plt.plot(frame_times, retroicor_regressors[:,0], label="sin1")
plt.plot(frame_times, 3 + retroicor_regressors[:,1], label="cos1")
plt.plot(frame_times, 6 + retroicor_regressors[:,2], label="sin2")
plt.plot(frame_times, 9 + retroicor_regressors[:,3], label="cos2")
plt.xlim([0,100])
plt.xlabel("time (s)", size=25)
plt.legend(ncol=2, prop={'size':12})
plt.tick_params(labelsize=20)
plt.yticks([])
pass<jupyter_output><empty_output><jupyter_text>Let's build now a design matrix that contains these regressors in addition to the motion parameters and an intercept.<jupyter_code>dm_retroicor = np.column_stack((np.ones(len(frame_times)), # Intercept
                                conf_df.filter(regex="rot|trans").to_numpy(), # Motion parameters
                                retroicor_regressors # ECG regressors
                               ))
dm_retroicor = pd.DataFrame(dm_retroicor, 
                            columns=["intercept", 
                                     "trans_x","trans_y","trans_z","rot_x","rot_y","rot_z",
                                     "sin1", "cos1", "sin2", "cos2"
                                      ]
                     )
dm_retroicor.index=frame_times
dm_retroicor.head()
fig, ax = plt.subplots()
ax.pcolormesh(dm_retroicor)
ax.set_xticks(0.5 + np.arange(dm_retroicor.shape[1]))
ax.set_xticklabels(dm_retroicor.columns, rotation=45, size=15)
pass<jupyter_output><empty_output><jupyter_text>Now fit this design matrix to our BOLD image<jupyter_code>first_level = FirstLevelModel(t_r=t_r, drift_model=None, signal_scaling=False, smoothing_fwhm=6)
first_level.fit(run_imgs=run_img, design_matrices=dm_retroicor)
retroicor_map = first_level.compute_contrast("sin1+cos1+sin2+cos2", stat_type="F", output_type="z_score")<jupyter_output><empty_output><jupyter_text>And plot this map in a glass brain fashion<jupyter_code>plotting.plot_glass_brain(retroicor_map, colorbar=True, threshold=1.96)<jupyter_output><empty_output><jupyter_text>Variation model's regressors Let's continue now with the generation of regressors using the Variation Model. The class for generating these regressors is found as `niphlem.models.HVPhysio`.<jupyter_code>from niphlem.models import HVPhysio<jupyter_output><empty_output><jupyter_text>Let's have a look at the documentation of the class<jupyter_code>print(HVPhysio.__doc__)
variaton_ecg = HVPhysio(physio_rate=400, t_r=t_r, delta=200, columns="mean", time_window=4.5,
                        low_pass=2.0, high_pass=0.5, peak_rise=0.5
                       )
variation_regressors = variaton_ecg.compute_regressors(signal=ecg_signal, time_scan=frame_times)
fig, ax = plt.subplots(figsize=(7,5))
plt.plot(frame_times, variation_regressors[:,0], label="ecg_var")
plt.legend()
plt.xlim([0,100])
plt.xlabel("time (s)", size=25)
plt.legend(prop={'size':12})
plt.tick_params(labelsize=20)
pass
dm_variation = np.column_stack((np.ones(len(frame_times)), # Intercept
                                conf_df.filter(regex="rot|trans").to_numpy(), # Motion parameters
                                variation_regressors # ECG regressors
                         ))

dm_variation = pd.DataFrame(dm_variation, 
                            columns=["intercept", 
                                     "trans_x","trans_y","trans_z","rot_x","rot_y","rot_z",
                                     "ecg_var"
                                    ]
                     )
dm_variation.index=frame_times
dm_variation.head()
fig, ax = plt.subplots()
ax.pcolormesh(dm_variation)
ax.set_xticks(0.5 + np.arange(dm_variation.shape[1]))
ax.set_xticklabels(dm_variation.columns, rotation=45, size=15)
pass<jupyter_output><empty_output><jupyter_text>Again, let's fit this design matrix to our BOLD data and show the pattern obtained<jupyter_code>first_level = FirstLevelModel(t_r=t_r, drift_model=None, signal_scaling=False, smoothing_fwhm=6)
first_level.fit(run_imgs=run_img, design_matrices=dm_variation)
variation_map = first_level.compute_contrast("ecg_var", stat_type="F", output_type="z_score")
plotting.plot_glass_brain(variation_map, colorbar=True, threshold=1.96)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09230326847503115, 'fraction_numerical': 0.023483178376305952, 'mean_word_length': 4.043015949734171, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 68, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '618098', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3555, 'n_tokens_neox': 3321, 'n_words': 899, 'path': 'examples/tutorial3.ipynb', 'repo_name': 'CoAxLab/brainhack-physio-project', 'size': '10433'}