This notebook walks through the preprocessing/QA workflow described in the course slides.
Importing modules
In [1]:
import os, errno
try:
datadir=os.environ['FMRIDATADIR']
assert not datadir==''
except:
datadir='/Users/poldrack/data_unsynced/myconnectome/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
from compute_fd_dvars import compute_fd,compute_dvars
from nipype.caching import Memory
mem = Memory(base_dir='.')
results_dir = os.path.abspath("../results")
if not os.path.exists(results_dir):
os.mkdir(results_dir)
def force_symlink(file1, file2):
try:
os.symlink(file1, file2)
except OSError, 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]:
We will need the mean EPI image later during the class so we will save it.
In [4]:
force_symlink(mcflirt_results.outputs.mean_img, os.path.join(results_dir, "meanbold.nii.gz"))
Use FSL's BET to obtain the brain mask
In [5]:
bet = mem.cache(fsl.BET)
bet_results = bet(functional=True,
in_file=mcflirt_results.outputs.mean_img,
mask=True)
bet_results.outputs
Out[5]:
Display the mean image, and show the outline of the brain mask.
In [6]:
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 [7]:
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
In [8]:
motiondata=numpy.loadtxt(mcflirt_results.outputs.par_file)
fd=compute_fd(motiondata)
dvars=compute_dvars(globalmean)
In [9]:
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 [10]:
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
Out[10]:
Check out the MELODIC report
In [11]:
from IPython.display import FileLink
FileLink(os.path.join(melodic_results.outputs.report_dir, '00index.html').split(os.getcwd() + '/')[1])
Out[11]:
Load the ICA components and examine their correlation with the motion signals. For components with a correlation > 0.4, show the component voxels.
In [12]:
ica_comps=numpy.loadtxt(os.path.join(melodic_results.outputs.out_dir,'melodic_mix'))
In [13]:
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]))
Take the first magnitue image and skull strip it
In [14]:
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 [15]:
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
Out[15]:
Prepare the fieldmap
In [16]:
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
Out[16]:
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)
Out[ ]:
Note that we don't really apply these transformations yet. We'll do it after first level modelling
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)
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)
Out[ ]:
High-pass filtering - this will take a few minutes
In [ ]:
hpfilt = mem.cache(fsl.maths.TemporalFilter)
hpcutoff = 120
TR = 1.16
hpfilt_results = hpfilt(highpass_sigma = hpcutoff/TR,
in_file=smooth_results.outputs.smoothed_file)
binarymaths = mem.cache(fsl.maths.BinaryMaths)
rescale_results = binarymaths(operation = 'add',
in_file=hpfilt_results.outputs.out_file,
operand_file=mcflirt_results.outputs.mean_img
)
Exercise: Plot the power spectrum for a voxel to look at the effects of filtering