Hi, I’m trying to fit some linear model to my data and get a contrast map using FitLins, and in the app’s arguments I see --space which seems to accept MNI space by default (line 140 here). I however couldn’t find any surface-based options such as fsaverage for this argument. Does the app support surface data too (can it read from fMRIPrep’s GIfTI files)? Thanks for the help!
Hi @sajjad,
It can currently work on CIFTI data: Usage — FitLins documentation. fsaverage is not supported yet. If you would like a Nilearn script that can work with GIFTIs (such as fsaverage or fsnative) let me know.
Best,
Steven
Thanks, @Steven! I was actually just going to move to Nilearn since it seemed like the second best option after FitLins, so yes I’d appreciate it if you could please share the script.
Cheers,
Sajjad
import os.path as op
import nibabel as nib
import nilearn.surface
from nilearn.glm.first_level import make_first_level_design_matrix, first_level_from_bids
from nilearn.glm.first_level.first_level import run_glm
from nilearn.glm.contrasts import compute_contrast
import numpy as np
import pandas as pd
from scipy import stats
from glob import glob
# Define parameters here
bids="/om2/user/smeisler/nsd_bids/" # Path to BIDS root
fmriprep_dir = op.join("derivatives","fmriprep") # BIDS-relative path to fMRIPrep
subject = "sub-01" # Subject name
task = "floc" # Task name
space = "fsnative" # BOLD projected on subject's freesurfer surface
hemi = "L" # L for left, R for right
use_smoothed = False
run_nums = ["1", "2", "3", "4", "5", "6"]
slice_time_ref = 0 # From the fMRIPrep command
# Final output dictionary for contrast objects (to be combined later)
contrast_objs = {}
# Loop over runs
for run_num in run_nums:
print("Processing run", run_num)
run = "_run-"+run_num # Run string in filename (define as empty string "" if no run label)
# Load GIFTI data and z-score it
#func_name = f"{subject}_task-{task}{run}_hemi-{hemi}_space-{space}_bold.func.gii"
func_name = f"{subject}_task-{task}{run}_space-{space}_hemi-{hemi}_bold.func.gii"
# If you smoothed data beforehand, make sure to point this to your smoothed file name!
if use_smoothed:
func_name = func_name.replace("_bold", "_desc-smoothed_bold")
gii_path = op.join(bids, fmriprep_dir, subject, "func", func_name)
gii_data = nilearn.surface.load_surf_data(gii_path)
gii_data_std = stats.zscore(gii_data,axis=1)
n_vertices= np.shape(gii_data)[0]
n_scans = np.shape(gii_data)[1]
# Use the volumetric data just to get the events and confounds file
img_filters = [('desc', 'preproc')]
if len(run)>0: img_filters.append(('run', run_num))
l1 = first_level_from_bids(
bids, task, space_label="T1w", sub_labels=[subject[4:]],
slice_time_ref=slice_time_ref,
hrf_model='spm',
drift_model=None, high_pass=None, drift_order=0, # Since we use fMRIPrep's cosine regressors
img_filters=img_filters,
memory_level=1, standardize=False, signal_scaling=0,
noise_model='ar1', verbose=0, n_jobs=1, minimize_memory=True,
derivatives_folder=fmriprep_dir)
# Extract information from the prepared model
t_r = l1[0][0].t_r
events = l1[2][0][0]
confounds = l1[3][0][0]
## From the confounds file, extract only those of interest
# Start with the motion and acompcor regressors
motion_keys = ['framewise_displacement', 'rot_x','rot_y','rot_z','trans_x','trans_y','trans_z']
wm_csf_keys = ['white_matter','csf']
# Now add non-steady-state volumes
non_steady_state_keys = [key for key in confounds.keys() if 'non_steady' in key]
# Add cosine regressors which act to high-pass filter data at 1/128 Hz
cosine_keys = [key for key in confounds.keys() if 'cosine' in key]
# Pull out the confounds we want to keep
confound_keys_keep = motion_keys + wm_csf_keys + cosine_keys + non_steady_state_keys
confounds_keep = confounds[confound_keys_keep]
# Set first value of FD column to the column mean
confounds_keep['framewise_displacement'][0] = np.nanmean(confounds_keep['framewise_displacement'])
## Add HRF models for our experimental conditions
# Start by getting times of scans
frame_times = t_r * (np.arange(n_scans) + slice_time_ref)
# Now use Nilearn to create the design matrix
design_matrix = make_first_level_design_matrix(frame_times,
events=events,
hrf_model='spm',
drift_model=None, # we use fMRIPrep's cosine regressors
add_regs=confounds_keep
)
# z-score the design matrix to standardize it
design_matrix_std = stats.zscore(design_matrix,axis=0)
# add constant back in to standardized design matrix since you cannot z-score a constant
design_matrix_std['constant'] = np.ones(len(design_matrix_std)).astype(int)
## Run the GLM
Y = np.transpose(gii_data_std)
X = np.asarray(design_matrix)
labels, estimates = run_glm(Y, X, verbose=1, n_jobs=-1)
## Define the contrast
contrast_matrix = np.eye(design_matrix.shape[1])
basic_contrasts = dict([(column, contrast_matrix[i])
for i, column in enumerate(design_matrix.columns)])
contrasts = {
'facesGTother': (
basic_contrasts['adult']/2
+ basic_contrasts['child']/2
- basic_contrasts['body']/8
- basic_contrasts['limb']/8
- basic_contrasts['number']/8
- basic_contrasts['word']/8
- basic_contrasts['car']/8
- basic_contrasts['instrument']/8
- basic_contrasts['corridor']/8
- basic_contrasts['house']/8
),
'wordsGTother': (
basic_contrasts['number']/2
+ basic_contrasts['word']/2
- basic_contrasts['body']/8
- basic_contrasts['limb']/8
- basic_contrasts['adult']/8
- basic_contrasts['child']/8
- basic_contrasts['car']/8
- basic_contrasts['instrument']/8
- basic_contrasts['corridor']/8
- basic_contrasts['house']/8
)
}
## Compute the contrast
for index, (contrast_id, contrast_val) in enumerate(contrasts.items()):
# Add a label to the output dictionary
if contrast_id not in contrast_objs:
contrast_objs[contrast_id] = []
# Define a template name for output statistical maps
outname_base = f"{subject}_task-{task}_contrast-{contrast_id}{run}_hemi-{hemi}_space-{space}_stat-X_statmap.func.gii" # generic name for stat
# compute contrast-related statistics
contrast = compute_contrast(labels, estimates, contrast_val,
contrast_type='t')
# add contrast to the output dictionary
contrast_objs[contrast_id].append(contrast)
# do the run-specific processing
betas = contrast.effect_size()
z_score = contrast.z_score()
t_value = contrast.stat()
p_value = contrast.p_value()
# Save the value maps as GIFTIs
# Effect size
gii_to_save = nib.gifti.gifti.GiftiImage()
gii_to_save.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=betas, datatype="NIFTI_TYPE_FLOAT32"))
outname = outname_base.replace("stat-X", "stat-effect")
if use_smoothed:
outname = outname.replace("_statmap", "_desc-smoothed_statmap")
nib.save(gii_to_save, outname)
# z-score
gii_to_save = nib.gifti.gifti.GiftiImage()
gii_to_save.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=z_score, datatype="NIFTI_TYPE_FLOAT32"))
outname = outname_base.replace("stat-X", "stat-z")
if use_smoothed:
outname = outname.replace("_statmap", "_desc-smoothed_statmap")
nib.save(gii_to_save, outname)
# t-value
gii_to_save = nib.gifti.gifti.GiftiImage()
gii_to_save.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=t_value, datatype="NIFTI_TYPE_FLOAT32"))
outname = outname_base.replace("stat-X", "stat-t")
if use_smoothed:
outname = outname.replace("_statmap", "_desc-smoothed_statmap")
nib.save(gii_to_save, outname)
# p-value
gii_to_save = nib.gifti.gifti.GiftiImage()
gii_to_save.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=p_value, datatype="NIFTI_TYPE_FLOAT32"))
outname = outname_base.replace("stat-X", "stat-p")
if use_smoothed:
outname = outname.replace("_statmap", "_desc-smoothed_statmap")
nib.save(gii_to_save, outname)
print("Producing Session-Wide Statistical Maps")
for index, (contrast_id, contrast_val) in enumerate(contrasts.items()):
contrast_concat = contrast_objs[contrast_id][0].__add__(contrast_objs[contrast_id][1]).__add__(contrast_objs[contrast_id][2]).__add__(contrast_objs[contrast_id][3]).__add__(contrast_objs[contrast_id][4]).__add__(contrast_objs[contrast_id][5])
betas = contrast_concat.effect_size()
z_score = contrast_concat.z_score()
t_value = contrast_concat.stat()
p_value = contrast_concat.p_value()
outname_base = f"{subject}_task-{task}_contrast-{contrast_id}_hemi-{hemi}_space-{space}_stat-X_statmap.func.gii" # generic name for stat
gii_to_save = nib.gifti.gifti.GiftiImage()
gii_to_save.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=betas, datatype="NIFTI_TYPE_FLOAT32"))
outname = outname_base.replace("stat-X", "stat-effect")
if use_smoothed:
outname = outname.replace("_statmap", "_desc-smoothed_statmap")
nib.save(gii_to_save, outname)
# z-score
gii_to_save = nib.gifti.gifti.GiftiImage()
gii_to_save.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=z_score, datatype="NIFTI_TYPE_FLOAT32"))
outname = outname_base.replace("stat-X", "stat-z")
if use_smoothed:
outname = outname.replace("_statmap", "_desc-smoothed_statmap")
nib.save(gii_to_save, outname)
# t-value
gii_to_save = nib.gifti.gifti.GiftiImage()
gii_to_save.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=t_value, datatype="NIFTI_TYPE_FLOAT32"))
outname = outname_base.replace("stat-X", "stat-t")
if use_smoothed:
outname = outname.replace("_statmap", "_desc-smoothed_statmap")
nib.save(gii_to_save, outname)
# p-value
gii_to_save = nib.gifti.gifti.GiftiImage()
gii_to_save.add_gifti_data_array(nib.gifti.gifti.GiftiDataArray(data=p_value, datatype="NIFTI_TYPE_FLOAT32"))
outname = outname_base.replace("stat-X", "stat-p")
if use_smoothed:
outname = outname.replace("_statmap", "_desc-smoothed_statmap")
nib.save(gii_to_save, outname)
1 Like
Note:
- that this loops over a set of runs and then makes a session-average map.
- I believe the
spaceandhemilabel order switched at some point (I am working with 20.2.7 outputs) so you may have to change that naming convention for thefunc_namevariable. - This uses volumetric output names to get the events and confounds files (since the gii files have weird naming conventions Nilearn is not set to handle natively through their BIDS grabber functions). I had
T1wspace images and used that for thespace_labelargument in thel1variable. If you have MNI you can put that there too. - Stats are done on one hemi at a time. If you want to produce both hemis you’ll have to run twice or add another layer to the loop.
- This code also is for a dataset without sessions. If you have session folders, you’ll have to update paths accordingly.
Hopefully the rest of it is straight forward and clear from the comments. If not, let me know!
Best,
Steven
Great, thanks so much! I’ll try this!
P.S. If I want to move towards group analysis, I think the best way would be through Nilearn’s glm.second_level.SecondLevelModel - @Steven have you tried this? Or, is there any other approach you’d recommend?
Best,
Sajjad
You can try running something similar to this (this was for the first level model). Y would be dimensioned something like num_vertices X n_subs, where the vertex values are the effect sizes from the first level models. X would be your second level design matrix of dimension n_subs X n_regressors.
Also, you’ll notice that that the first level results are not smoothed for this case, but if you want to smooth your giis, you’ll want to make sure those files are being loaded by the script.
Here is an example smoothing script with FreeSurfer:
#!/bin/bash -l
bids="/path/to/nsd_bids/" # Or replace with your own BIDS dataset
subject="sub-01" # Or replace with your own subject ID
task="floc" # Which task to process
run="1" # Which run to process
hemi="L" # Hemi name in fMRIPrep outputs; "L" or "R"
if [ "$hemi" == "L" ];
then hemi_fs="lh"; # Hemi name in FreeSurfer conventions; "lh" or "rh"
elif [ "$hemi" == "R" ];
then hemi_fs="rh"; # Hemi name in FreeSurfer conventions; "lh" or "rh"
fi;
# Hemi name in FreeSurfer conventions; "lh" or "rh"
space="fsnative" # Using native space surface outputs
fwhm="4" # Desired smoothing kernel size (mm FWHM)
fmriprep_dir=${bids}/derivatives/fmriprep/
freesurfer_dir=${bids}/derivatives/freesurfer/
export SUBJECTS_DIR=${freesurfer_dir} # Tell FreeSurfer where subjects live
gii_in="${fmriprep_dir}/${subject}/func/${subject}_task-${task}_run-${run}_space-${space}_hemi-${hemi}_bold.func.gii" # Input file name
gii_out=${gii_in/_bold/_desc-smoothed_bold} # Add smoothed label for output name
mris_fwhm --i ${gii_in} --o ${gii_out} --so \
--fwhm ${fwhm} --subject ${subject} --hemi ${hemi_fs} # Perform the smoothing
1 Like