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)
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')
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
Out[3]:
In [4]:
fsl.MCFLIRT.help()
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
Out[6]:
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)
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()
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)
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.
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 - 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 [ ]: