Quality control task activation analysis
The interactive version of this notebook is available at https://github.com/TheAxonLab/hcph-sops/tree/mkdocs/docs/analysis
Load data and extract parameters¶
import nibabel as nb
import numpy as np
import pandas as pd
from pathlib import Path
import matplotlib.pyplot as plt
Locate all the QCT functional scans and corresponding files¶
Our dataset, called the Human Connectome Phantom (HCPh) dataset, contains 36 sessions and each session has been acquired on a single healthy individual with the same acquisition protocol. One of the scans acquired during those sessions is a BOLD task fMRI. The task named Quality-Control Task (QCT) comprises four paradigms: a central fixation dot (blank), gaze movement, visual grating patterns, and finger-tapping (left/right). We ran our dataset through fMRIPrep to preprocess the data and extract confounds that we will use to denoise the data. We also generated the events file that lists the onset and duration of each stimuli from the task files output by Psychopy. Let us locate all the files that we need using the convenient BIDS interface implemented in Nilearn.
# Nilearn works well with file paths and there's a nice BIDS interface implemented
from nilearn.interfaces import bids as nibids
# This will work if the "fmriprep" dir is a valid BIDS dataset
# (requires to have "sub-###" directories and a "dataset_description.json")
path_to_dataset = Path("/data/datasets/hcph")
path_to_derivatives = Path("/data/derivatives/hcph-derivatives/fmriprep-23.2.2-wp")
# Load the T1w images so we can plot activity on top
t1w = path_to_derivatives / "sub-001/anat/sub-001_acq-undistorted_desc-preproc_T1w.nii.gz"
# Load the QCT functional scans
qct_files = nibids.get_bids_files(
path_to_derivatives, file_tag="bold", file_type="nii.gz", filters=[("task", "qct")]
)
#Keep only the files with no space entity, that is the file in BOLD native space (filtering with ("space", None) did not work on the previous line)
qct_files = [file for file in qct_files if 'space-' not in file]
# Load the QCT events
qct_events_all = nibids.get_bids_files(
path_to_dataset, file_tag="events", file_type="tsv", filters=[("task", "qct")]
)
# Filter out pilot and excluded sessions, as well as session for which fMRIPrep did not run
qct_events = []
for element in qct_events_all:
if (
"ses-0" in element
):
qct_events.append(element)
"""
# Code useful to understand why we do not have the same number of events than files
import re
match_events = set([re.search(r'ses-0\d\d', s).group() if re.search(r'ses-0\d\d', s) else None for s in qct_events])
match_files = set([re.search(r'ses-0\d\d', s).group() if re.search(r'ses-0\d\d', s) else None for s in qct_files])
print(match_events.symmetric_difference(match_files))
"""
print(len(qct_events))
print(len(qct_files))
36 36
Extract acquisition parameters¶
For building the model of the BOLD during task, we need to extract the timing of acquisition which depends on the TR and the number of volumes, so let us extract those parameters using nibabel.
# Extract acquisition parameters
tr = []
n_scans = []
for file in qct_files:
fmri = nb.load(file)
tr.append(fmri.header.get_zooms()[3])
n_scans.append(fmri.shape[3])
Quality control of acquisition parameters¶
In our single-site dataset, all sessions have been acquired at the same site with the same acquisition protocol. As such, we have to check that all functional scans have the same repetition time (TR) and the same number of volumes. A deviation in TR or number of volumes has to be investigated further to understand what happened and the scan might need to be excluded.
print(tr)
print(n_scans)
[1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6, 1.6] [99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99]
Good news! All the scans have the same TR and the same number of volumes. Let's simplify notation by storing only one tr
and one n_samples
value. Additionally, now that those two parameters are defined, we can extract the acquisition timing.
tr = tr[0]
n_scans = n_scans[0]
frame_times = np.arange(n_scans) * tr
Extract confounds¶
fMRI data are notoriously noisy, as such we need to perform denoising through a process known as nuisance regression. fMRIPrep conveniently compute a series of regressors that you might want to remove from your signal. Thus, let us load some confounds from fMRIPrep derivatives for using them later.
# Nilearn also has a nice fMRIPrep interface implemented that we will use
# to load some confounds extracted by the tool
from nilearn.interfaces import fmriprep
qct_confounds, qct_sample_masks = fmriprep.load_confounds(
qct_files,
strategy=["motion", "wm_csf", "scrub"],
motion="basic",
wm_csf="basic",
scrub=5,
fd_threshold=0.4,
std_dvars_threshold=1.5,
)
Extract activation maps for each subject¶
Put together design matrix¶
To get information about what regions of the brain are activated during a certain task, you must first know when the task of interest was happening during the scan. The information about the timing of each task is stored as columns (one per task plus other derivatives) in a matrix called design matrix. We use nilearn to compute the design matrix from the events files.
from nilearn.glm.first_level import make_first_level_design_matrix
hrf_model = "glover"
design_matrices = []
for filename, confounds in zip(qct_events, qct_confounds):
# Load the event file as a dataframe
event = pd.read_csv(filename, sep="\t")
# Remove eye-tracker related event as they are of no interest here
event = event[(event['trial_type'] != 'et-record-on') & (event['trial_type'] != 'et-record-off')]
# Construct design matrix
design_matrices.append(
make_first_level_design_matrix(
frame_times, event, hrf_model=hrf_model, add_regs=confounds
)
)
/home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn( /home/cprovins/.local/lib/python3.11/site-packages/nilearn/glm/first_level/experimental_paradigm.py:167: UserWarning: The following unexpected columns in events data will be ignored: value warnings.warn(
Note the warning that has been emitted by the make_first_level_design_matrix
: "UserWarning: The following unexpected columns in events data will be ignored: value", we will come back to it later.
Let us now have a look at one design matrix as an example. In a event-related design like here, the design matrices are constructed by convolving a Haemodynamic Response Function (HRF) with timeseries indicative of tasks' onset and duration. As such, we expect peaks around the time of stimuli presentation with a smooth transition around that peak.
from nilearn.plotting import plot_design_matrix
plt.figure(figsize=(3, 3))
plot_design_matrix(design_matrices[0])
plt.title("ses-1 design matrix", fontsize=12)
# Adjust layout to prevent clipping of titles and labels
plt.tight_layout()
# Show the plot
plt.show()
<Figure size 300x300 with 0 Axes>
Note that we indicated the model to take into account confounds (e.g. motion parameters rot*
and trans*
), so they appear in the design matrix alongside the task regressors (blank
,cog
, mot
and vis
).
It seems that the task regressors do present the pattern we were refering to, but it is not very clear. A clearer way to visualize this is to plot the task regressors as timeseries.
import re
# Quick check of the event models in signal plots
plt.figure(figsize=(3, 3))
fig, axs = plt.subplots(12, 3, figsize=(24, 70))
events_sessions = [
re.search(r"ses-0\d\d", s).group() if re.search(r"ses-0\d\d", s) else None
for s in qct_events
]
# Plot task regressors for each QCT scans
for i, (design_matrix, ses) in enumerate(zip(design_matrices, events_sessions)):
space_between_lines = 2
# Just plot the task regressors not the confounds
dm_subset = design_matrix[["blank", "cog", "mot", "vis"]]
for j, column in enumerate(dm_subset):
axs.flat[i].plot(
dm_subset.loc[:, column] + j * space_between_lines, label=column
)
axs.flat[i].set_title(ses, fontsize=18)
axs.flat[i].set_yticks(np.arange(dm_subset.shape[1]) * space_between_lines)
axs.flat[i].grid()
axs.flat[i].legend(fontsize=18)
# Adjust layout to prevent clipping of titles and labels
plt.tight_layout()
# Show the plot
ax = plt.show()
<Figure size 300x300 with 0 Axes>