This notebook walks through the preprocessing/QA workflow described in the course slides.

Importing modules


In [1]:
import os, errno, sys
try:
    datadir=os.environ['FMRIDATADIR']
    assert not datadir==''
except:
    datadir='/home/vagrant/data/ds031/sub00001'
print('Using data from',datadir)

%matplotlib inline

from nipype.interfaces import fsl, nipy
import nibabel
import numpy
import nilearn.plotting
import matplotlib.pyplot as plt
sys.path.append('/home/vagrant/fmri-analysis-vm/analysis/utils')
from compute_fd_dvars import compute_fd,compute_dvars
from nipype.caching import Memory
mem = Memory(base_dir='.')
results_dir = os.path.abspath("/home/vagrant/fmri-analysis-vm/results")
if not os.path.exists(results_dir):
    os.mkdir(results_dir)
    
def force_symlink(file1, file2):
    try:
        os.symlink(file1, file2)
    except OSError as e:
        if e.errno == errno.EEXIST:
            os.remove(file2)
            os.symlink(file1, file2)


Using data from /home/vagrant/data

Setting up variables


In [2]:
subject='ses014'  

bolddir=os.path.join(datadir,'ds031/sub00001',subject,
        'functional')
boldfile=os.path.join(bolddir,'sub00001_ses014_task002_run001_bold.nii.gz')

Motion correction using MCFLIRT

This will take about a minute.


In [3]:
mcflirt = mem.cache(fsl.MCFLIRT)
mcflirt_results = mcflirt(in_file=boldfile,
                          save_plots=True,
                          mean_vol=True)
mcflirt_results.outputs


161006-17:31:13,7 workflow INFO:
	 Executing node be24d7d11ebf3374aa93d29f4bbedd75 in dir: /home/vagrant/fmri-analysis-vm/analysis/preprocessing/nipype_mem/nipype-interfaces-fsl-preprocess-MCFLIRT/be24d7d11ebf3374aa93d29f4bbedd75
161006-17:31:13,11 workflow INFO:
	 Collecting precomputed outputs
Out[3]:
mat_file = <undefined>
mean_img = /home/vagrant/fmri-analysis-vm/analysis/preprocessing/nipype_mem/nipype-interfaces-fsl-preprocess-MCFLIRT/be24d7d11ebf3374aa93d29f4bbedd75/sub00001_ses014_task002_run001_bold_mcf.nii.gz_mean_reg.nii.gz
out_file = /home/vagrant/fmri-analysis-vm/analysis/preprocessing/nipype_mem/nipype-interfaces-fsl-preprocess-MCFLIRT/be24d7d11ebf3374aa93d29f4bbedd75/sub00001_ses014_task002_run001_bold_mcf.nii.gz
par_file = /home/vagrant/fmri-analysis-vm/analysis/preprocessing/nipype_mem/nipype-interfaces-fsl-preprocess-MCFLIRT/be24d7d11ebf3374aa93d29f4bbedd75/sub00001_ses014_task002_run001_bold_mcf.nii.gz.par
rms_files = <undefined>
std_img = <undefined>
variance_img = <undefined>

In [4]:
fsl.MCFLIRT.help()


Wraps command **mcflirt**

Use FSL MCFLIRT to do within-modality motion correction.

For complete details, see the `MCFLIRT Documentation.
<http://www.fmrib.ox.ac.uk/fsl/mcflirt/index.html>`_

Examples
--------
>>> from nipype.interfaces import fsl
>>> from nipype.testing import example_data
>>> mcflt = fsl.MCFLIRT(in_file=example_data('functional.nii'), cost='mutualinfo')
>>> res = mcflt.run() # doctest: +SKIP

Inputs::

	[Mandatory]
	in_file: (an existing file name)
		timeseries to motion-correct
		flag: -in %s, position: 0

	[Optional]
	args: (a unicode string)
		Additional parameters to the command
		flag: %s
	bins: (an integer (int or long))
		number of histogram bins
		flag: -bins %d
	cost: ('mutualinfo' or 'woods' or 'corratio' or 'normcorr' or
		 'normmi' or 'leastsquares')
		cost function to optimize
		flag: -cost %s
	dof: (an integer (int or long))
		degrees of freedom for the transformation
		flag: -dof %d
	environ: (a dictionary with keys which are a bytes or None or a value
		 of class 'str' and with values which are a bytes or None or a value
		 of class 'str', nipype default value: {})
		Environment variables
	ignore_exception: (a boolean, nipype default value: False)
		Print an error message instead of throwing an exception in case the
		interface fails to run
	init: (an existing file name)
		inital transformation matrix
		flag: -init %s
	interpolation: ('spline' or 'nn' or 'sinc')
		interpolation method for transformation
		flag: -%s_final
	mean_vol: (a boolean)
		register to mean volume
		flag: -meanvol
	out_file: (a file name)
		file to write
		flag: -out %s
	output_type: ('NIFTI_PAIR_GZ' or 'NIFTI' or 'NIFTI_PAIR' or
		 'NIFTI_GZ')
		FSL output type
	ref_file: (an existing file name)
		target image for motion correction
		flag: -reffile %s
	ref_vol: (an integer (int or long))
		volume to align frames to
		flag: -refvol %d
	rotation: (an integer (int or long))
		scaling factor for rotation tolerances
		flag: -rotation %d
	save_mats: (a boolean)
		save transformation matrices
		flag: -mats
	save_plots: (a boolean)
		save transformation parameters
		flag: -plots
	save_rms: (a boolean)
		save rms displacement parameters
		flag: -rmsabs -rmsrel
	scaling: (a float)
		scaling factor to use
		flag: -scaling %.2f
	smooth: (a float)
		smoothing factor for the cost function
		flag: -smooth %.2f
	stages: (an integer (int or long))
		stages (if 4, perform final search with sinc interpolation
		flag: -stages %d
	stats_imgs: (a boolean)
		produce variance and std. dev. images
		flag: -stats
	terminal_output: ('stream' or 'allatonce' or 'file' or 'none')
		Control terminal output: `stream` - displays to terminal immediately
		(default), `allatonce` - waits till command is finished to display
		output, `file` - writes output to file, `none` - output is ignored
	use_contour: (a boolean)
		run search on contour images
		flag: -edge
	use_gradient: (a boolean)
		run search on gradient images
		flag: -gdt

Outputs::

	mat_file: (a list of items which are an existing file name)
		transformation matrices
	mean_img: (an existing file name)
		mean timeseries image
	out_file: (an existing file name)
		motion-corrected timeseries
	par_file: (an existing file name)
		text-file with motion parameters
	rms_files: (a list of items which are an existing file name)
		absolute and relative displacement parameters
	std_img: (an existing file name)
		standard deviation image
	variance_img: (an existing file name)
		variance image

References::

We will need the mean EPI image later during the class so we will save it.


In [5]:
force_symlink(mcflirt_results.outputs.mean_img, os.path.join(results_dir, "meanbold.nii.gz"))
force_symlink(mcflirt_results.outputs.par_file, os.path.join(results_dir, "motion.par"))

Use FSL's BET to obtain the brain mask


In [6]:
bet = mem.cache(fsl.BET)
bet_results = bet(functional=True,
              in_file=mcflirt_results.outputs.mean_img,
              mask=True)
force_symlink(bet_results.outputs.mask_file, os.path.join(results_dir, "mask.nii.gz"))
bet_results.outputs


161006-17:31:13,54 workflow INFO:
	 Executing node 91d85fb6ab4be0b07acdfb9b82d6f6ac in dir: /home/vagrant/fmri-analysis-vm/analysis/preprocessing/nipype_mem/nipype-interfaces-fsl-preprocess-BET/91d85fb6ab4be0b07acdfb9b82d6f6ac
161006-17:31:13,58 workflow INFO:
	 Collecting precomputed outputs
Out[6]:
inskull_mask_file = <undefined>
inskull_mesh_file = <undefined>
mask_file = /home/vagrant/fmri-analysis-vm/analysis/preprocessing/nipype_mem/nipype-interfaces-fsl-preprocess-BET/91d85fb6ab4be0b07acdfb9b82d6f6ac/sub00001_ses014_task002_run001_bold_mcf.nii.gz_mean_reg_brain_mask.nii.gz
meshfile = <undefined>
out_file = /home/vagrant/fmri-analysis-vm/analysis/preprocessing/nipype_mem/nipype-interfaces-fsl-preprocess-BET/91d85fb6ab4be0b07acdfb9b82d6f6ac/sub00001_ses014_task002_run001_bold_mcf.nii.gz_mean_reg_brain.nii.gz
outline_file = <undefined>
outskin_mask_file = <undefined>
outskin_mesh_file = <undefined>
outskull_mask_file = <undefined>
outskull_mesh_file = <undefined>
skull_mask_file = <undefined>

Display the mean image, and show the outline of the brain mask.


In [7]:
mask_display=nilearn.plotting.plot_epi(mcflirt_results.outputs.mean_img,cmap='gray')
mask_display.add_contours(bet_results.outputs.mask_file, levels=[.5])


Compute and plot the global signal within the mask across timepoints


In [ ]:
maskdata=nibabel.load(bet_results.outputs.mask_file).get_data()
bolddata=nibabel.load(mcflirt_results.outputs.out_file).get_data()
maskvox=numpy.where(maskdata)
globalmean=numpy.zeros(bolddata.shape[3])
globalcv=numpy.zeros(bolddata.shape[3])
for t in range(bolddata.shape[3]):
    tmp=bolddata[:,:,:,t]
    globalmean[t]=numpy.mean(tmp[maskvox])
    globalcv[t]=numpy.std(tmp[maskvox])/numpy.mean(tmp[maskvox])

del bolddata

Load motion data and compute FD/DVARS.

Exercise: The next cell will open the motion-corrected fMRI file in fslview. Explore the data at the timepoints near the maximum point of motion (which is output by the next cell). What do you see in this image and those adjacent to it? Did motion correction fully remove the effects of motion in those timepoints?


In [ ]:
motiondata=numpy.loadtxt(mcflirt_results.outputs.par_file)

fd=compute_fd(motiondata)
dvars=compute_dvars(globalmean)
print('Maximum motion at timepoint %d'%numpy.where(fd==numpy.max(fd))[0])

_=os.system('fslview %s'%mcflirt_results.outputs.out_file)


Maximum motion at timepoint 107

In [ ]:
f, (ax1, ax2,ax3,ax4)=plt.subplots(4, sharex=True)
ax1.plot(globalmean)
ax1.set_title('Mean global signal for in-mask voxels')

ax2.plot(globalcv)
ax2.set_title('Coefficient of variation for in-mask voxels')

ax3.plot(fd)
ax3.set_title('Framewise Displacement')

ax4.plot(dvars)
ax4.set_title('DVARS')
plt.tight_layout()

Performing Independent Component Analysis to look for artefacts

Run independent components analysis on the data using MELODIC - this will take a few minutes.


In [ ]:
melodic = mem.cache(fsl.MELODIC)
melodic_results = melodic(out_all=True,
                          report=True,
                          in_files=mcflirt_results.outputs.out_file,
                          mask=bet_results.outputs.mask_file)
melodic_results.outputs

Check out the MELODIC report


In [ ]:
from IPython.display import FileLink
FileLink(os.path.join(melodic_results.outputs.report_dir, '00index.html').split(os.getcwd() + '/')[1])

Load the ICA components and examine their correlation with the motion signals. For components with a correlation > 0.4, show the component voxels.


In [ ]:
ica_comps=numpy.loadtxt(os.path.join(melodic_results.outputs.out_dir,'melodic_mix'))
os.listdir(melodic_results.outputs.out_dir)

In [ ]:
ica_motion_corr=numpy.zeros(ica_comps.shape[1])
for c in range(ica_comps.shape[1]):
    ica_motion_corr[c]=numpy.corrcoef(ica_comps[:,c],fd)[0,1]
    if abs(ica_motion_corr[c])>0.4:
        comp_img=nibabel.load(os.path.join(melodic_results.outputs.out_dir,'stats/thresh_zstat%d.nii.gz'%int(c+1)))
        nilearn.plotting.plot_stat_map(comp_img,mcflirt_results.outputs.mean_img,threshold=1.5,
               title='component %d - r(FD)=%0.3f'%(c+1,ica_motion_corr[c]))

In [ ]:
import seaborn as sns
sns.distplot(ica_motion_corr)

Fieldmap unwarping

Take the first magnitue image and skull strip it


In [ ]:
magfile=os.path.join(datadir,'ds031/sub00001',subject,
                'fieldmap/sub00001_ses014_001_magnitude.nii.gz')
roi = mem.cache(fsl.ExtractROI)
pick_first_mag = roi(in_file=magfile, t_min=0, t_size=1)

# we'll need this later
force_symlink(pick_first_mag.outputs.roi_file, os.path.join(results_dir, "fmapmag.nii.gz"))

mag_bet_results = bet(functional=True,
                      in_file=pick_first_mag.outputs.roi_file,
                      mask=True)
mag_bet_results.outputs
mask_display=nilearn.plotting.plot_epi(pick_first_mag.outputs.roi_file,cmap='gray')
mask_display.add_contours(mag_bet_results.outputs.mask_file, levels=[.5])

This was not quite tight enough. Let's try with other parameters


In [ ]:
mag_bet_results = bet(functional=True,
                      in_file=pick_first_mag.outputs.roi_file,
                      mask=True,
                      frac=0.65)
force_symlink(mag_bet_results.outputs.out_file, os.path.join(results_dir, "fmapmagbrain.nii.gz"))
mask_display=nilearn.plotting.plot_epi(pick_first_mag.outputs.roi_file,cmap='gray')
mask_display.add_contours(mag_bet_results.outputs.mask_file, levels=[.5])
mag_bet_results.outputs

Prepare the fieldmap


In [ ]:
prepare = mem.cache(fsl.PrepareFieldmap)
prepare_results = prepare(in_phase = os.path.join(datadir,'ds031/sub00001',subject,
                                                  'fieldmap/sub00001_ses014_001_phasediff.nii.gz'),
                          in_magnitude = mag_bet_results.outputs.out_file,
                          output_type = "NIFTI_GZ")
force_symlink(prepare_results.outputs.out_fieldmap, os.path.join(results_dir, "fieldmap.nii.gz"))
prepare_results.outputs

Which parts of the brain were most deformed by B0 inhomogeneities?


In [ ]:
nilearn.plotting.plot_stat_map(nilearn.image.index_img(prepare_results.outputs.out_fieldmap,0), 
                               pick_first_mag.outputs.roi_file, threshold=100)

Note that we don't really apply these transformations yet. We'll do it after first level modelling


In [ ]:
fugue = mem.cache(fsl.FUGUE)
fugue_results = fugue(in_file = mcflirt_results.outputs.mean_img,
                      fmap_in_file = prepare_results.outputs.out_fieldmap,
                      unwarp_direction = "y",
                      dwell_time = 2.6/10000.0,
                      mask_file = mag_bet_results.outputs.mask_file)
fugue_results.outputs

Exercise: Compute the difference between the original mean EPI and the field unwarped one, and display the difference map.

Spatial smoothing

Perform spatial smoothing using Gaussian kernel. Can you spot the difference between top and bottom rows?


In [ ]:
smooth = mem.cache(fsl.utils.Smooth)
smooth_results = smooth(fwhm=2.5,
                        in_file=mcflirt_results.outputs.out_file,
                        output_type = "NIFTI")

nilearn.plotting.plot_epi(mcflirt_results.outputs.mean_img,colorbar=True, vmin=410, vmax=21000)
nilearn.plotting.plot_epi(nilearn.image.mean_img(smooth_results.outputs.smoothed_file), colorbar=True, vmin=410, vmax=21000)

High pass filtering

High-pass filtering - this will take a few minutes


In [ ]:
hpfilt = mem.cache(fsl.maths.TemporalFilter)
TR = 2.32
hpfilt_results = hpfilt(highpass_sigma = 100/(2*TR),
                        in_file=smooth_results.outputs.smoothed_file,
                        output_type = "NIFTI")

mean = mem.cache(fsl.maths.MeanImage)
mean_results = mean(in_file = smooth_results.outputs.smoothed_file)

rescale = mem.cache(fsl.maths.BinaryMaths)
rescale_results = rescale(in_file=hpfilt_results.outputs.out_file,
                          operand_file = mean_results.outputs.out_file,
                          operation = "add",
                          output_type = "NIFTI")

force_symlink(rescale_results.outputs.out_file, os.path.join(results_dir, "preprocessed_epi_native_space.nii"))

Exercise: Plot the times series for a voxel before and after filtering. What is the effect?


In [ ]: