NSB Pipeline for fMRI

NSB Pipeline for fMRI

fMRI can be tricky to preprocess, depending on the information available and the fMRI acquisition sequence. We have put together some basic recommendations below:

<style> </style>

SAMPLE SIZE	SINGLE BAND	SINGLE ECHO	PIPELINE
Small (e.g. <100 subjects)	Yes	Yes	fMRIPrep + ICA-AROMA + Friston24R + WMSR +CSFSR +GMSR + Physiological Noise Filtering + bandpass filtering
		No	fMRIPrep + tedana + Friston24R + WMSR +CSFSR +GMSR + Physiological Noise Filtering + bandpass filtering
	No	Yes	fMRIPrep (STC off) + ICA-AROMA + Friston24R + WMSR +CSFSR +GMSR + Physiological Noise Filtering + bandpass filtering
		No	fMRIPrep (STC off) + tedana + Friston24R + WMSR +CSFSR +GMSR+ Physiological Noise Filtering + bandpass filtering
Large (e.g. >100 subjects)	Yes	Yes	fMRIPrep + ICA-FIX + Friston24R + WMSR +CSFSR +GMSR + Physiological Noise Filtering + bandpass filtering
		No	fMRIPrep + tedana + Friston24R + WMSR +CSFSR +GMSR + Physiological Noise Filtering + bandpass filtering
	No	Yes	fMRIPrep (STC off) + ICA-FIX + Friston24R + WMSR +CSFSR +GMSR + Physiological Noise Filtering + bandpass filtering
		No	fMRIPrep (STC off) + tedana + Friston24R + WMSR +CSFSR +GMSR+ Physiological Noise Filtering + bandpass filtering

The very first step would be to make sure your data is in BIDS format for fMRIPrep to recognise. BIDSification + fMRIPrep is always the first step. If you have multiband data, turn slice timing off in fMRIPrep. If your data is multi-echo, you will need to raise the --me-output-echos flag for independent echo preprocessing.

SPECIAL NOTE ON AROMA/FIX AFTER PREP: fMRIPrep runs mcflirt for us already (yay!) - and this is required by AROMA/FIX. However, it is hidden in the work folder. We can point to this file when running AROMA, or copy it into the MELODIC/FEAT directory prior to running FIX.

e.g.

/home/tconstab1/kg98_scratch/Toby/derivatives/fMRIPrep/rest/work/fmriprep_23_1_wf/single_subject_${s}_wf/func_preproc_task_rest_dir_RL_echo_1_wf/bold_hmc_wf/mcflirt/sub-${s}_*mcf.nii.gz.par

CAUTION

If you are running FIX after fMRIPrep, MAKE SURE ALL OF THE DENOISING STRATEGIES ARE TURNED OFF IN MELODIC/FEAT. This is incredibly important. Doubling up on data preprocessing can lead to some very odd/strange/plain wrong consequences. Also, if you smooth by hand yourself prior to MELODIC/FEAT, make sure the spatial smoothing is set to 0.

SPECIAL NOTE ON SPATIAL SMOOTHING: fMRIPrep does not apply any spatial smoothing in the workflow it implements. tedana developers do not recommend spatial smoothing, so you will not really need to worry about adding any intermediary steps in a multi-echo pipeline between the fMRIPrep and tedana workflows (see more on using tedana after fMRIPrep in the tedana section of the wiki). Special care needs to be paid if applying ICA-FIX or ICA-AROMA, though. Spatial smoothing to a degree of ~3mm-6mm is generally recommended prior to either of these approaches. Historically, we have generally estimated noise components on the spatially smoothed data, and then stripped the noise components using fsl_regfilt:

fsl_regfilt -i original_nonsmoothed_data -o denoised_data -d filtered_func_data.ica/melodic_mix -f "2,5,9"

One can then apply the WM/CSF/GM SRs at the same time as stripping noise ICs (better for preserving degrees of freedom) on the NON smoothed data, or after stripping the noise ICs from the NON smoothed fMRI img (may be better for removal of any residual noise not captured during ICA denoising). The latter approach tends to be more widely practiced.

To perform WM / CSF/ GM signal regression using fMRIPrep outputs, follow the following steps:

1. Resample the masks to BOLD space (first volume) using the transforms available.

output_dir=/home/tconstab1/kg98_scratch/Toby/WHOLEMBBP/workspace/PipelineComparisons/ResampleT1TisseuMaskstoNative
freesurfer_dir=/home/tconstab1/kg98_scratch/Toby/WHOLEMBBP/workspace/derivatives/fmriPREP_pepolar/rest

cd $output_dir

if [ ! -d ${output_dir}/sub-${s} ]; then
	mkdir sub-${s}
	fslroi ${freesurfer_dir}/sub-${s}/func/sub-${s}_task-rest_dir-RL_echo-2_desc-preproc_bold.nii.gz sub-${s}/sub-${s}_firstvolume_task-rest_dir-RL_echo-2_desc-preproc_bold.nii.gz 0 1
	for TISSUE in CSF WM GM
	do
	antsApplyTransforms -d 3 -i ${freesurfer_dir}/sub-${s}/anat/sub-${s}_label-${TISSUE}_probseg.nii.gz -r $output_dir/sub-${s}/sub-${s}_firstvolume_task-rest_dir-RL_echo-2_desc-preproc_bold.nii.gz -o sub-${s}/sub-${s}_label-${TISSUE}_probseg_Native.nii.gz -n Linear -t ${freesurfer_dir}/sub-${s}/func/sub-${s}_task-rest_dir-RL_from-T1w_to-scanner_mode-image_xfm.txt
	done
fi

Also check that the transforms worked by overlay using fsleyes.

2. Erode the CSF / WM masks to avoid partial voluming effects. We tend to erode the masks up to 5 times - selecting the latest erosion iteration that maintains at least 5 voxels in the mask. We do not worry about this for the GM mask.

    		#Erode CSF mask once
    		fslmaths ${csf_probability_map} -ero -kernel 3D sub-${s}_space-NATIVE_label-CSF_probseg_eroded_1_times.nii.gz
    		#Erode wm mask 5x
    		input_mask=${wm_probability_map}
    		for ((i = 1; i < 6; i++)); do
        			output_mask="sub-${s}_space-NATIVE_label-WM_probseg_eroded_${i}_times.nii.gz"
        			fslmaths "${input_mask}" -ero -kernel 3D "${output_mask}"
        			input_mask="${output_mask}"
    		done
			# Find the file with most erosions and at least 5 voxels
			selected_wm_mask=""
			for ((i = 5; i >= 1; i--)); do
    				mask="sub-${s}_space-NATIVE_label-WM_probseg_eroded_${i}_times.nii.gz"
    				voxel_count=$(fslstats "${mask}" -V | awk '{print $1}')
            
    				if [ "${voxel_count}" -ge 5 ]; then
        				selected_wm_mask="${mask}"
        				break
    				fi
			done

			# Define filenames for the best eroded WM mask and eroded CSF mask
			mv ${selected_wm_mask} best_erosion_${selected_wm_mask}

			#The CSF erosion sometimes deletes ALL voxels. If so, use original mask
			selected_eroded_csf_mask=sub-${s}_space-NATIVE_label-CSF_probseg_eroded_1_times.nii.gz
            voxel_count_CSF=$(fslstats "sub-${s}_space-NATIVE_label-CSF_probseg_eroded_1_times.nii.gz" -V | awk '{print $1}')
			selected_eroded_csf_mask=${selected_eroded_csf_mask}
    		if [ "${voxel_count_CSF}" -ge 5 ]; then
        		selected_eroded_csf_mask=${csf_probability_map}
    		fi
			mv ${selected_eroded_csf_mask} best_erosion_${selected_eroded_csf_mask}

3. Extract the timeseries from the NON smoothed ICA denoised image using fslmeants

fslmeants -i $input_file_aroma -m ${csf_probability_map} -o ${confounds_output}/sub-${s}_average_csf_signal.txt
fslmeants -i $input_file_aroma -m ${wm_probability_map} -o ${confounds_output}/sub-${s}_average_wm_signal.txt
fslmeants -i $input_file_aroma -m ${gm_probability_map} -o ${confounds_output}/sub-${s}_average_gm_signal.txt

4. Find the differences between successive points (dx[i] = x[i] - x[i-1]) to approximate the numerical derivatives of these tissue specific timeseries vectors (same method as in fMRIPrep to derive tissue specific timeseries vector derivatives - but they do not erode their masks)

5. Calculate the squares of the raw tissue specific timeseries vectors and their derivatives

6. Merge the Friston 24 motion parameter vectors (3 translation params, 3 rotation params, each with derivatives, powers, and power derivatives) from the fMRIPrep Confounds .tsv file with the raw, derivative and power derivative WM/CSF/GM timeseries vectors. Perform regression on the NON smoothed ICA denoised image e.g. using fsl_regfilt (NOTE - other methods of regression are available - e.g. clean_img and signal_clean in python - these tend to result in images that look nothing like an fMRI image! Myself and others prefer fsl for this reason).

NOW, the final step is to filter out frequencies associated with non-neuronal signal.

to do this you do a bandpass filter between the frequencies you want to filter out. This can be done with the following script with the following inputs.

matlab -nodisplay -r "bandpass_filter('$functional_file_loc', '"$functional_file', '$brain_mask', '$TR', '$lp_hz', '$hp_hz'); exit;"

lp_hz is often set to 0.008, while hp_hz is set to 0.08. note that it requires the image processing toolbox and Rest_Idealfilter

function bandpass_filter(func_folder, file_name, brain_mask, TR, lp, hp)
% This function bandpass filters a 4D niftiimage, needs a brain mask 3D
% file so that only brain voxels are filtered.
TR=str2double(TR);
lp=str2double(lp);
hp=str2double(hp);

to_load=append(func_folder,'/',file_name);
data=niftiread(to_load);
mask=niftiread(brain_mask);

% Keep only brain voxels for filtering
mask=reshape(mask,[],1);
mask=logical(mask);

dim=size(data);
data=reshape(data,[],dim(4));
data=data';
numVoxels=size(data,2);
time=size(data,1);
data=data(:,mask);

[filtered_dat] = rest_IdealFilter(data,TR,[lp hp]);

% Put filtered brain back into full image and reshape back to 4D
data_out=zeros(time,numVoxels);
data_out(:,mask)=filtered_dat;
data_out=data_out';
data_out=reshape(data_out,dim);
data_out=single(data_out);

name_out=extractBefore(file_name,".");
name_out=[name_out,'_bpass.nii'];
file_out=string([func_folder,'/',name_out]);
niftiwrite(data_out,file_out)
gzip([func_folder,'/',name_out]);
delete([func_folder,'/',name_out]);
end

Congratulations ! You're all done.

The last step is fMRI QC. Please refer to the fMRI QC wiki page for more details.