Evaluation of fMRIPrep: comparison with FSL feat

This notebook is a supplemental material to the paper: [doi here]

0. Setting up

First, make sure the FMRIPREP_DATA_HOME environment variable is set and pointing to the root folder of the BIDS structure for OpenfMRI's DS000030.

In [1]:
/home/oesteban/.anaconda3/lib/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters
In [2]:
def mean_std_map(pipe_home, meanmask, force=False, lazy=False, maskval=1000):
    pipe_std = pipe_home / 'summary_stdev.nii.gz'
    pipe_mean = pipe_home / 'summary_means.nii.gz'

    if force or not pipe_mean.is_file():
        all_mus = []
        if lazy:
            all_mus = [nb.load(str(f)) for f in pipe_home.glob(
        if not all_mus:
            print('Generating means file')
            pipe_files = list(pipe_home.glob(
            all_mus = []
            for f in pipe_files:
                mean = mean_img(str(f))
                data = mean.get_data()
                sigma = np.percentile(data[meanmask], 50) / maskval
                data /= sigma
                all_mus.append(nb.Nifti1Image(data, mean.affine, mean.header))
        meannii = concat_imgs(all_mus, auto_resample=False)
        force = True

    if force or not pipe_std.is_file():
        print('Generating standard deviation map')
        meannii = nb.load(str(pipe_mean))
        nb.Nifti1Image(meannii.get_data().std(3), meannii.affine, meannii.header).to_filename(str(pipe_std))
    return pipe_mean, pipe_std

# Use the WM mask to normalize intensities of EPI means
meanmask = nb.load(str(ATLAS_HOME / '2mm_tpm_wm.nii.gz')).get_data() > 0.9

1. Standard deviation across subjects of each subject's BOLD average.

We normalize the BOLD averages to have the same median value within the WM mask (WM mask of the MNI template)

In [3]:
# Calculate average and std
fprep_mean, fprep_std = mean_std_map(fprep_home, meanmask)
feat_mean, feat_std = mean_std_map(feat_home, meanmask)

coords = [-27, 0, 7]

thres = 20
vmin = 50
vmax = 200

# Plot
fig = plt.figure(figsize=(20,10))

plotting.plot_anat(str(feat_std), cut_coords=coords, colorbar=True, cmap='cividis',
                   threshold=thres, vmin=vmin, vmax=vmax, title='feat',
plotting.plot_anat(str(fprep_std), cut_coords=coords, colorbar=True, cmap='cividis',
                   threshold=thres, vmin=vmin, vmax=vmax, title='fmriprep',
plotting.plot_glass_brain(str(feat_std), threshold=200, colorbar=True, title='feat',
plotting.plot_glass_brain(str(fprep_std), threshold=200, colorbar=True, title='fmriprep',
<Figure size 432x288 with 0 Axes>

2. First-level analysis

2.1. Inspecting one subject

In [4]:
pipelines = ['fslfeat','fmriprep']
subject = 'sub-11090'
home = {'fslfeat': feat_home, 'fmriprep': fprep_home}

Let's have a look to the z-statistic map for both pipelines

In [5]:
images = {}
bg_images = {}
preptpl = '{0}_task-stopsignal_bold_space-MNI152NLin2009cAsym_preproc.nii.gz'.format
for pipeline in pipelines:
    z11 = ANALYSIS_HOME / subject / 'func' / '{}_task-stopsignal_variant-{}_zstat11.nii.gz'.format(subject, pipeline)
    images[pipeline] = str(z11)
    im = nb.load(str(home[pipeline] / subject / 'func' / preptpl(subject)))
    bg_images[pipeline] = nb.Nifti1Image(im.get_data().mean(3), im.affine, im.header)

# Plot
fig = plt.figure(figsize=(20,10))

i = 1
for idx, pipeline in enumerate(pipelines):
    plotting.plot_stat_map(images[pipeline], title=pipeline, cmap='RdYlBu_r',
                           colorbar=True,symmetric_cbar=True, alpha=0.6,
                           axes=plt.subplot(2, 2, i))
    i += 1
                              axes=plt.subplot(2, 2, i))
    i += 1

<Figure size 432x288 with 0 Axes>

Here we plot the distribution of z-values for this particular subject

In [6]:

fprep_vals = nb.load(images['fmriprep']).get_data()[mask2mm]
feat_vals = nb.load(images['fslfeat']).get_data()[mask2mm]
sn.distplot(fprep_vals[np.abs(fprep_vals) > 1.68], label='fmriprep', kde=False, norm_hist=True)
sn.distplot(feat_vals[np.abs(feat_vals) > 1.68], label='feat', kde=False, norm_hist=True)
plt.title("Distribution of Z-values of 1st level analysis - GO-StopSuccess contrast")
<matplotlib.legend.Legend at 0x7f11e099b828>

These z-values distributions correspond to the following activation maps (for this one particular subject)

In [7]:
cut_coords = [-15, -8, 6, 30, 46, 62]
for idx,pipeline in enumerate(pipelines):
   plotting.plot_stat_map(images[pipeline], title=pipeline, vmax=5,
                          display_mode='z', threshold=1.65, cut_coords=cut_coords,