# - 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 common modules
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 appending the directory containing the module to the sys.path list.
Import the sys module, and show the path variable from the sys module.
# Show path in sys module
import sys
sys.path
Find the directory where you put the spm_funcs.py module. For example, for me, that would be /Users/mb312/pna_code. Append this directory to the path list.
# Append the directory containing `spm_funcs.py` to sys.path
sys.path.append('/Users/mb312/pna_code')
# Check you can import spm_funcs
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 file has one line of text 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
task = np.loadtxt('ds114_sub009_t2r1_cond.txt')
task.shape
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
ons_durs = task[:, :2] / 2.5
ons_durs
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')
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 = np.zeros(img.shape[-1])
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. That means the task started at the beginning of scan index 4, and lasted for 12 scans. That means there should be 12 consecutive 1 values in time_course, starting at index 4. That means that index 4 + 12 = 16 should be zero again, because there are 12 values starting at (including 4) going up to (but not including) 16.
# - try running this if you don't believe me
len(range(4, 16))
So, for the first row, you will want to set time_course[4] through time_course[15] equal to 1.
# Fill in values of 1 for positions of on blocks in time course
for onset, duration in ons_durs:
time_course[onset:onset + duration] = 1
Plot the time course:
# Plot the time course
plt.plot(time_course)
Comparing task to rest¶
Make a boolean array is_task_tr which is True when time_course is 1 and False otherwise.
Make another array is_rest_tr that is the opposite - True when time_course is 0 and False otherwise.
# Make two boolean arrays encoding task, rest
is_task_tr = (time_course == 1)
is_rest_tr = (time_course == 0)
Now read the image data into an array:
# Read the image data into an array.
data = img.get_data()
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
on_volumes = data[..., is_task_tr]
Now select the volumes where the time course is 0 (rest volumes):
# Create a new 4D array only containing the rest volumes
off_volumes = data[..., is_rest_tr]
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.
on_mean = on_volumes.mean(axis=-1)
off_mean = off_volumes.mean(axis=-1)
Now subtract the rest mean from the task mean to get a difference volume.
# Create a difference volume
difference = on_mean - off_mean
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(difference[:, :, 14], 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 the outlier easier to see
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.subplots to show a matching slice from the bad volume and the following volume side-by-side:
# Select bad volume and following volume
# Use subplots to show slice from good and next
bad_vol = data[..., 0]
next_vol = data[..., 1]
fig, axes = plt.subplots(1, 2)
axes[0].imshow(bad_vol[:, :, 14], cmap='gray')
axes[1].imshow(next_vol[:, :, 14], cmap='gray')
In fact, 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
off_volumes_fixed = off_volumes[..., 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
off_mean_fixed = off_volumes_fixed.mean(axis=-1)
difference_fixed = on_mean - off_mean_fixed
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
plt.subplot(1, 2, 1)
plt.imshow(difference[:, :, 14], cmap='gray')
plt.subplot(1, 2, 2)
plt.imshow(difference_fixed[:, :, 14], cmap='gray')
If you are racing ahead¶
We can actually do this calculation of the mean image in another way, by doing an elementwise array multiplication, followed by a sum.
See if you can work out how to do that. This section is difficult, so don't be surprised if you have to stare at it for a long time, or you can't get it easily. The exercise is to practice slicing and thinking about 4D arrays.
Imagine that we are looking at a single voxel time series, for 3D voxel coordinate $i, j, k$.
Call the time series vector $y_0, y_1 ... y_{T-1}$ where
- $y_0 = image\_data[i, j, k, 0]$,
- $y_1 = image\_data[i, j, k, 1]$
- $...$
- $y_{T-1} = image\_data[i, j, k, T-1]$
$T = 173$ - the number of TRs.
Call the indices of the rest volumes sequence $R \in 0, 1, 2, 3, 16 ... $ of length $N$, and the indices of the task volumes $K \in 4, 5, 6 ...$ of length $M$.
The average of the points from $y$ within the rest volumes is:
$$ \frac{1}{N} \sum_r{y_r} : r \in R $$The average of the points from $y$ within the task volumes is:
$$ \frac{1}{M} \sum_k{y_k} : k \in K $$Subtracting rest from task to get the difference in means at this voxel:
$$ \frac{1}{M} \sum_k{y_k} - \frac{1}{N} \sum_r{y_r} $$Moving the division inside the sums:
$$ \sum_k{\frac{1}{M} y_k} - \sum_r{\frac{1}{N} y_r} $$We can join these two sums by iterating over the entire time series $y_0, y_1 ... y_{T-1}$, selecting the weight to apply according to whether the index $t \in 0, 1, ... T-1$ is a rest scan index or a task index:
$$ \sum_t{ S y_t } \begin{cases} S = -\frac{1}{N} ,& \text{if } t \in R \\ S = \frac{1}{M} ,& \text{if } t \in K \\ \end{cases} $$- Make a new 4D array for
Sin the equation above, the same shape as the image data array (hint: usenp.zeros); - Set the values for the array
Sto $-1/N$ for all the values in the rest volumes and $1/M$ for all the values in the task volumes (hint: use your boolean slicing arrays); - Do the $Sy_t$ calculation for every voxel and time point (hint: elementwise multiplication);
- Do the sum (hint:
axisagain); - Check you get the same answer as you did when you subtracted the means before.
S = np.zeros(data.shape)
N = np.sum(is_rest_tr)
M = np.sum(is_task_tr)
S[..., is_rest_tr] = -1. / N
S[..., is_task_tr] = 1. / M
new_difference = np.sum(data * S, axis=-1)
plt.imshow(new_difference[:, :, 14], cmap='gray')
For extra points, can you change the S array so that S removes the bad volume from the difference calculation?
S[..., is_rest_tr] = -1. / (N - 1)
S[..., 0] = 0
new_fixed_difference = np.sum(data * S, axis=-1)
plt.imshow(new_fixed_difference[:, :, 14], cmap='gray')
