# Compatibility with Python 3
from __future__ import print_function # print('me') instead of print 'me'
from __future__ import division # 1/2 == 0.5, not 0
# Show figures inside the notebook
%matplotlib inline
import numpy as np # the Python array package
import matplotlib.pyplot as plt # the Python plotting package
Using our module¶
Later on we are going to use our new module spm_funcs.py.
In order to do this, we need to make sure Python can find the module. Python stores all the directories that it will search for modules in the path variable in the sys module. We can make sure Python can find our new module by inserting the directory containing the module into sys.path.
Import the sys module, and show the path variable from the sys module.
# Show path in sys module
Find the directory where you put the spm_funcs.py module. For example, for me, that would be /Users/mb312/pna_code. Use the insert method of the path variable to insert this directory at the top of the path. Check the help for the insert method for how to do this.
# Insert the directory containing `spm_funcs.py` at the top
# of sys.path
# Check you can import spm_funcs
import sys
sys.path.insert(0,"/home/jb/aaa/pna_code")
import spm_funcs
The task time-course¶
We've previously had a look at the file ds114_sub009_t2r1.nii. This is a 4D FMRI image.
Now we want to see whether we can detect any signal in that image relating to the task.
The task was a block design with 10 seconds rest followed by 7 repeats of (30 seconds when the subject thought of verbs followed by 30 seconds rest).
There's a file in this directory called ds114_sub009_t2r1_cond.txt with the task information in it. It comes from subject 9 task 2 and run 1 from OpenFMRI dataset https://openfmri.org/dataset/ds000114 .
The has one line per "on" block, giving the onset time (in seconds), the duration (in seconds) and the amplitude (expected amount of activation for this block).
Here are the first four lines:
10 30.000000 1
70 30.000000 1
130 30.000000 1
190 30.000000 1# Read the file into an array called "task".
# "task" should have 3 columns (onset, duration, amplitude)
# HINT: np.loadtxt
tasktime = np.loadtxt("ds114_sub009_t2r1_cond.txt")
print(tasktime)
The TR for this FMRI run was 2.5 seconds. We need to convert the onsets and durations to TRs - so for example the first onset was at 10 seconds, which was the start of TR 4 (10 / 2.5).
Select out the first two columns of task, and divide by the TR to convert the onset and duration times to be in terms of TRs:
# Select first two columns and divide by TR
tasktime[:,:2] /= 2.5
print(tasktime)
Our next step is to make an on-off vector that is 0 when the subject is doing nothing and 1 when the subject is doing the covert verb task.
The vector will have one value (either 0 or 1) for each TR.
First use nibabel to load the image ds114_sub009_t2r1.nii. Check the image shape to find the number of TRs.
# Load the image and check the image shape to get the number of TRs
import nibabel as nib
img = nib.load('ds114_sub009_t2r1.nii')
print(img.shape)
Next make a new vector called time_course with one entry per TR, with all elements in the vector being zero:
# Make new zero vector
time_course = zeros((img.shape[-1]))
print(time_course.shape)
Loop over the rows in the onsets / durations array to give you an onset / duration pair. For each of these pairs, set the matching positions in the time_course vector to 1. For example, the first pair will be 4, 12. For that row, you will want to set time_course[4] through time_course[15] equal to 1.
# Fill in 1s for positions of on blocks in time course
for t0, dur in tasktime[:,:2]:
time_course[round(t0):round(t0+dur)] = 1.
print(time_course)
Plot the time course:
# Plot the time course
plt.plot(time_course)
Comparing task to rest¶
Make a boolean array which is True when time_course is 1 and False when time_course is 0.
Make another array that is the opposite - False when time_course is 0 and True when time course is 1.
# Make two boolean arrays encoding task, rest
time_course = time_course.astype(bool)
opposite_time_course = (~time_course)
plt.plot(opposite_time_course)
Now read the image data into an array:
# Read the image data into an array.
imgdata = img.get_data()
imgdata.shape
Remember that the 4D array consists of one volume (3D array) per TR.
We want to select the volumes where the time course is 1 (task volumes). Do this by slicing, using the boolean array you just made.
# Create a new 4D array only containing the task volumes
taskvol = imgdata[...,time_course]
taskvol.shape
Now select the volumes where the time course is 0 (rest volumes):
# Create a new 4D array only containing the rest volumes
restvol = imgdata[...,opposite_time_course]
restvol.shape
We want to know whether there is a difference in signal in the task volumes compared to the rest volumes. Take the mean over the task volumes and mean over the rest volumes. You should end up with two 3D volumes.
# Create the mean volume across all the task volumes
# Then create the mean volume across all the rest volumes
# Hint: remember the `axis` keyword.
Now subtract the rest mean from the task mean to get a difference volume.
# Create a difference volume
diffvol = taskvol.mean(axis=3) - restvol.mean(axis=3)
diffvol.shape
Show a slice over the third dimension of the difference volume, from somewhere around the center of the third axis:
# Show a slice over the third dimension
plt.imshow(diffvol[:,:,15],cmap='gray')
This is the difference between activation and rest. It looks a little strange. Maybe there are some artefacts here.
Fixing an artefact¶
Now let's use our spm_globals function to see if there are any volumes with odd values for this estimate of the mean.
Use the get_spm_globals function in spm_funcs to get the SPM global values for each volume in ds114_sub009_t2r1.nii.
# Get SPM global values for each volume
spm_globals = spm_funcs.get_spm_globals('ds114_sub009_t2r1.nii')
# Plot the SPM global values for each volume
# Try using 'o' marker symbol to make it easier to see, as in
# plt.plot(variable, 'o')
plt.plot(spm_globals,'o')
There is one volume that is a large outlier in the SPM globals.
Slice out the bad volume (to a 3D array) and the following volume.
Use plt.subplot to show a matching slice from the bad volume and the following volume side-by-side:
# Select bad volume and following volume
# Use subplot to show slice from good and next
f, ax = subplots(1,2)
ax[0].imshow(restvol[:,:,15,0])
ax[1].imshow(restvol[:,:,15,1])
As you might be able to see, the bad volume is one of the rest volumes.
Use your slicing skills to remove this volume from your selection of rest volumes.
# Use slicing to remove outlier volume from rest volumes
restvol_corrected = restvol[...,1:]
Make a new mean for the rest volumes, and subtract this mean from the mean for the task volumes to make a new difference image.
Give the new difference image a new name, so we can compare to the old difference image later.
# Make new mean for rest volumes, subtract from task mean
diffvol_corr = taskvol.mean(axis=3) - restvol_corrected.mean(axis=3)
Use subplot to show the same slice from the old difference volume and the new difference volume, side by side:
# Use subplot to show same slice from old and new difference volume
f, ax = subplots(1,2)
ax[0].imshow(diffvol[:,:,15])
ax[1].imshow(diffvol_corr[:,:,15])
