Hello! Thanks in advance for any advice/help!
Summary of what happened:
I’ve been running first and second-level GLMs using Nilearn to create task contrast maps for functional MRI n-back task data. I’ve tried multiple different strategies (including cosine regressors/aCompCor components, and using a less strict high-pass filter), and every time I run the second-level model I get a lot of deactivation in the 0back-baseline contrast, which is quite unexpected. The 2back-baseline and 2back-0back contrasts look relatively normal. The most recent strategy I used was a conservative high pass filter (0.005 Hz), regular signal scaling, and non_steady_state in the confounds strategy. Does anyone have any idea on how to troubleshoot/what might be causing these deactivations in the 0back-baseline contrast?
First-level GLM script:
#!/usr/bin/env python
"""Run first-level GLM for n-back task within a single subject."""
import sys
from pathlib import Path
import numpy as np
import pandas as pd
from scipy.stats import norm
from nilearn.glm.first_level import first_level_from_bids
from nilearn.interfaces.bids import save_glm_to_bids
# ============================================================
# GET SUBJECT FROM COMMAND LINE
# ============================================================
sub_id = sys.argv[1] # subject passed by SLURM array / CLI
sub_labels = [sub_id] # nilearn expects a list
print(f"Running first-level GLM for subject: {sub_id}")
# BIDS root
bids_root = Path("/cbica/projects/executive_function/data/bids/EF_bids_data_DataLad")
# fMRIPrep derivatives containing preprocessed func + confounds
derivatives_folder = (
"/cbica/projects/executive_function/EF_dataset/derivatives/"
"fmriprep_unzipped/fmriprep_func"
)
deriv_root = Path(derivatives_folder)
# Task and space
task_label = "nback"
space_label = "MNI152NLin6Asym"
# ============================================================
# FIND FMRIPREP BRAIN MASK FOR THIS SUBJECT
# ============================================================
# Example expected filename:
# sub-21050_ses-1_task-nback_run-01_space-MNI152NLin6Asym_res-2_desc-brain_mask.nii.gz
# or sub-20259_ses-3_task-nback_acq-VARIANTObliquity_run-01_space-MNI152NLin6Asym_res-2_desc-brain_mask.nii.gz
mask_pattern = (
f"sub-{sub_id}/ses-*/func/"
f"sub-{sub_id}_ses-*_task-{task_label}*"
f"_space-{space_label}_res-2_desc-brain_mask.nii.gz"
)
mask_candidates = sorted(deriv_root.glob(mask_pattern))
mask_img = None
if len(mask_candidates) == 0:
print(
f"[WARNING] No fMRIPrep brain mask found for sub-{sub_id} with pattern:\n"
f" {mask_pattern}\n"
"Falling back to Nilearn's automatic mask.\n"
)
else:
# XXX: One run per session, not subject, right?
# Since only one run per subject, take the first match
mask_img = str(mask_candidates[0])
print(f"Using fMRIPrep brain mask for sub-{sub_id}:\n {mask_img}\n")
# ============================================================
# BUILD FIRST-LEVEL MODEL FROM BIDS
# ============================================================
# This model applies across sessions within a subject.
MODEL_TYPES = ["rtdur", "nortdur"]
for model_type in MODEL_TYPES:
# Output directory for maps + reports
out_dir = Path(f"/cbica/projects/executive_function/code/task_contrast/final/results_troubleshoot4/first-level/nback-{model_type}")
out_dir.mkdir(parents=True, exist_ok=True)
(
models,
models_run_imgs,
models_events,
models_confounds,
) = first_level_from_bids(
dataset_path=bids_root,
task_label=task_label,
space_label=space_label,
sub_labels=sub_labels,
derivatives_folder=derivatives_folder,
img_filters=[("desc", "preproc")],
smoothing_fwhm=5.0,
n_jobs=4,
verbose=1,
# Use fMRIPrep's cosine regressors for drift (high_pass),
# so we turn OFF the GLM's own drift model.
drift_model="cosine",
high_pass=0.005,
# Confound strategy: motion + high_pass (cosines) + aCompCor (first 5)
confounds_strategy=("motion", "non_steady_state"),
confounds_motion="basic",
#confounds_compcor="anat_combined",
#confounds_n_compcor=5,
# Use fMRIPrep mask if found; otherwise None (Nilearn auto-mask)
mask_img=mask_img,
)
# Single subject
model = models[0]
run_imgs = models_run_imgs[0]
events_list = models_events[0]
confounds_list = models_confounds[0]
updated_events_list = []
for events in events_list:
events = events.copy()
# Rename trial_type labels before building design matrix
if "trial_type" in events.columns:
events["trial_type"] = events["trial_type"].replace(
{
"0BACK": "zero_back",
"2BACK": "two_back",
}
)
if model_type == "rtdur":
# Implement Jeanette Mumford's ConsDurRTDur model
# and rename trial types to valid Python identifiers:
# 0BACK -> zero_back
# 2BACK -> two_back
# Create dataframe only containing trials with responses
response_events = events.loc[~np.isnan(events["response_time"])].copy()
# Set duration to response time
response_events.loc[:, "duration"] = response_events.loc[:, "response_time"]
# Change trial type to a single value (RTDur)
response_events.loc[:, "trial_type"] = "RTDur"
# Add new "condition" back into events dataframe
events = pd.concat((events, response_events))
events = events.sort_values(by="onset")
updated_events_list.append(events)
print(f"Number of runs for subject {sub_id}: {len(run_imgs)}")
print("Confounds entries for each run:")
for i, c in enumerate(confounds_list):
print(f" Run {i}: {c}")
# ---------- FIT MODEL ----------
model.minimize_memory = False
print(f"\nFitting {model_type} GLM for subject {sub_id}...")
model = model.fit(
run_imgs,
events=updated_events_list,
confounds=confounds_list,
)
# Inspect design matrix
design_matrix = model.design_matrices_[0]
print("\nDesign matrix columns:")
print(design_matrix.columns)
print(f"\nTotal # regressors in design matrix: {design_matrix.shape[1]}")
save_glm_to_bids(
model,
contrasts=["two_back - zero_back", "two_back", "zero_back"],
contrast_types={
"two_back - zero_back": "t",
"two_back": "t",
"zero_back": "t",
},
out_dir=out_dir,
threshold=norm.isf(0.001),
cluster_threshold=10,
bg_img=(
"/cbica/projects/executive_function/templateflow/"
"tpl-MNI152NLin6Asym/tpl-MNI152NLin6Asym_res-02_T1w.nii.gz"
),
)
print("\nDone.\n")
Second-level GLM script:
#!/usr/bin/env python
"""Run second-level GLM for n-back task across all subjects."""
from pathlib import Path
import numpy as np
import pandas as pd
from nilearn.glm.second_level import SecondLevelModel
from nilearn.image import load_img
from nilearn.interfaces.bids import save_glm_to_bids
CONTRAST_LABELS = ["twoBack", "zeroBack", "twoBackMinusZeroBack"]
MODEL_TYPES = ["rtdur", "nortdur"]
group_mask_img = load_img(
"/cbica/projects/executive_function/.cache/templateflow/"
"tpl-MNI152NLin6Asym/tpl-MNI152NLin6Asym_res-02_desc-brain_mask.nii.gz"
)
bg_img = load_img(
"/cbica/projects/executive_function/.cache/templateflow/"
"tpl-MNI152NLin6Asym/tpl-MNI152NLin6Asym_res-02_desc-brain_T1w.nii.gz"
)
for model_type in MODEL_TYPES:
firstlevel_dir = Path(f"/cbica/projects/executive_function/code/task_contrast/final/results_troubleshoot4/first-level/nback-{model_type}")
for contrast_label in CONTRAST_LABELS:
print("\n====================================================")
print(f"Running group model for contrast: {contrast_label}")
print("====================================================\n")
group_out_dir = Path(
f"/cbica/projects/executive_function/code/task_contrast/final/results_troubleshoot4/second-level/nback-{model_type}/group-{contrast_label}"
)
group_out_dir.mkdir(exist_ok=True, parents=True)
print(f"MODEL TYPE: {model_type}")
print(f"Reading first-level maps from: {firstlevel_dir}")
print(f"Writing second-level outputs to: {group_out_dir}")
pattern = (
"sub-*/sub-*_task-nback_space-MNI152NLin6Asym_"
f"contrast-{contrast_label}_stat-effect_statmap.nii.gz"
)
effect_maps = sorted(firstlevel_dir.glob(pattern))
if len(effect_maps) == 0:
raise RuntimeError(
f"No first-level maps found for contrast '{contrast_label}' with pattern:\n {pattern}"
)
print(f"Found {len(effect_maps)} maps before QC.")
# ----- QC: drop maps that are all-nan or constant -----
good_maps = []
good_subs = []
for p in effect_maps:
img = load_img(p)
data = img.get_fdata()
if not np.isfinite(data).any():
print(f"[DROP] {p} is all non-finite.")
continue
if np.nanstd(data) == 0:
print(f"[DROP] {p} is constant/zero.")
continue
good_maps.append(p)
good_subs.append(p.name.split("_")[0].replace("sub-", ""))
effect_maps = good_maps
subject_labels = good_subs
print(f"Kept {len(effect_maps)} maps after QC.")
if len(effect_maps) < 2:
raise RuntimeError(
f"Need at least 2 valid subjects for group model of {contrast_label} "
f"(found {len(effect_maps)})."
)
# Design matrix
design_matrix = pd.DataFrame(
{"intercept": [1.0] * len(effect_maps)},
index=subject_labels,
)
# Fit second-level
model = SecondLevelModel(mask_img=group_mask_img, minimize_memory=False)
model = model.fit(effect_maps, design_matrix=design_matrix)
# Save results
contrasts = {contrast_label: "intercept"}
save_glm_to_bids(
model=model,
contrasts=contrasts,
out_dir=group_out_dir,
threshold=0.001,
height_control="fpr",
cluster_threshold=10,
bg_img=bg_img,
two_sided=True,
)
print(f"\nSaved second-level outputs to:\n {group_out_dir}\n")
Version:
nilearn version 0.12.1
Environment (Docker, Singularity / Apptainer, custom installation):
I am running the software on a high performance computing cluster in an environment where I have nilearn installed. I installed the software using pip install.
