Last but not least, the 2nd-level analysis. After we normalized all subject data into template space, we can now do the group analysis. To show the flexibility of Nipype, we will run the group analysis on data with two different smoothing kernel (fwhm= [4, 8]) and two different normalization (ANTs and SPM).
This example will also directly include thresholding of the output, as well as some visualization.
Let's start!
In [ ]:
%pylab inline
from os.path import join as opj
from nipype.interfaces.io import SelectFiles, DataSink
from nipype.interfaces.spm import (OneSampleTTestDesign, EstimateModel,
EstimateContrast, Threshold)
from nipype.interfaces.utility import IdentityInterface
from nipype.pipeline.engine import Workflow, Node
from nipype.interfaces.fsl import Info
from nipype.algorithms.misc import Gunzip
In [ ]:
experiment_dir = '/output'
output_dir = 'datasink'
working_dir = 'workingdir'
# Smoothing withds used during preprocessing
fwhm = [4, 8]
# Which contrasts to use for the 2nd-level analysis
contrast_list = ['con_0001', 'con_0002', 'con_0003', 'con_0004', 'con_0005',
'ess_0006', 'ess_0007']
# Mask to use for the group analysis
mask = Info.standard_image('MNI152_T1_2mm_brain_mask.nii.gz')
In [ ]:
# Gunzip - unzip the mask image
gunzip = Node(Gunzip(in_file=mask), name="gunzip")
# OneSampleTTestDesign - creates one sample T-Test Design
onesamplettestdes = Node(OneSampleTTestDesign(),
name="onesampttestdes")
# EstimateModel - estimates the model
level2estimate = Node(EstimateModel(estimation_method={'Classical': 1}),
name="level2estimate")
# EstimateContrast - estimates group contrast
level2conestimate = Node(EstimateContrast(group_contrast=True),
name="level2conestimate")
cont1 = ['Group', 'T', ['mean'], [1]]
level2conestimate.inputs.contrasts = [cont1]
# Threshold - thresholds contrasts
level2thresh = Node(Threshold(contrast_index=1,
use_topo_fdr=True,
use_fwe_correction=False,
extent_threshold=0,
height_threshold=0.005,
height_threshold_type='p-value',
extent_fdr_p_threshold=0.05),
name="level2thresh")
In [ ]:
# Infosource - a function free node to iterate over the list of subject names
infosource = Node(IdentityInterface(fields=['contrast_id', 'fwhm_id']),
name="infosource")
infosource.iterables = [('contrast_id', contrast_list),
('fwhm_id', fwhm)]
# SelectFiles - to grab the data (alternativ to DataGrabber)
templates = {'cons': opj(output_dir, 'norm_spm', 'sub-*_fwhm{fwhm_id}',
'w{contrast_id}.nii')}
selectfiles = Node(SelectFiles(templates,
base_directory=experiment_dir,
sort_filelist=True),
name="selectfiles")
# Datasink - creates output folder for important outputs
datasink = Node(DataSink(base_directory=experiment_dir,
container=output_dir),
name="datasink")
# Use the following DataSink output substitutions
substitutions = [('_contrast_id_', '')]
subjFolders = [('%s_fwhm_id_%s' % (con, f), 'spm_%s_fwhm%s' % (con, f))
for f in fwhm
for con in contrast_list]
substitutions.extend(subjFolders)
datasink.inputs.substitutions = substitutions
In [ ]:
# Initiation of the 2nd-level analysis workflow
l2analysis = Workflow(name='l2analysis')
l2analysis.base_dir = opj(experiment_dir, working_dir)
# Connect up the 2nd-level analysis components
l2analysis.connect([(infosource, selectfiles, [('contrast_id', 'contrast_id'),
('fwhm_id', 'fwhm_id')]),
(selectfiles, onesamplettestdes, [('cons', 'in_files')]),
(gunzip, onesamplettestdes, [('out_file',
'explicit_mask_file')]),
(onesamplettestdes, level2estimate, [('spm_mat_file',
'spm_mat_file')]),
(level2estimate, level2conestimate, [('spm_mat_file',
'spm_mat_file'),
('beta_images',
'beta_images'),
('residual_image',
'residual_image')]),
(level2conestimate, level2thresh, [('spm_mat_file',
'spm_mat_file'),
('spmT_images',
'stat_image'),
]),
(level2conestimate, datasink, [('spm_mat_file',
'2ndLevel.@spm_mat'),
('spmT_images',
'2ndLevel.@T'),
('con_images',
'2ndLevel.@con')]),
(level2thresh, datasink, [('thresholded_map',
'2ndLevel.@threshold')]),
])
In [ ]:
# Create 1st-level analysis output graph
l2analysis.write_graph(graph2use='colored', format='png', simple_form=True)
# Visualize the graph
from IPython.display import Image
Image(filename=opj(l2analysis.base_dir, 'l2analysis', 'graph.dot.png'))
Out[ ]:
In [ ]:
l2analysis.run('MultiProc', plugin_args={'n_procs': 4})
In [ ]:
# Change the SelectFiles template and recreate the node
templates = {'cons': opj(output_dir, 'norm_ants', 'sub-*_fwhm{fwhm_id}',
'{contrast_id}_trans.nii')}
selectfiles = Node(SelectFiles(templates,
base_directory=experiment_dir,
sort_filelist=True),
name="selectfiles")
# Change the substituion parameters for the datasink
substitutions = [('_contrast_id_', '')]
subjFolders = [('%s_fwhm_id_%s' % (con, f), 'ants_%s_fwhm%s' % (con, f))
for f in fwhm
for con in contrast_list]
substitutions.extend(subjFolders)
datasink.inputs.substitutions = substitutions
Now, we just have to recreate the workflow.
In [ ]:
# Initiation of the 2nd-level analysis workflow
l2analysis = Workflow(name='l2analysis')
l2analysis.base_dir = opj(experiment_dir, working_dir)
# Connect up the 2nd-level analysis components
l2analysis.connect([(infosource, selectfiles, [('contrast_id', 'contrast_id'),
('fwhm_id', 'fwhm_id')]),
(selectfiles, onesamplettestdes, [('cons', 'in_files')]),
(gunzip, onesamplettestdes, [('out_file',
'explicit_mask_file')]),
(onesamplettestdes, level2estimate, [('spm_mat_file',
'spm_mat_file')]),
(level2estimate, level2conestimate, [('spm_mat_file',
'spm_mat_file'),
('beta_images',
'beta_images'),
('residual_image',
'residual_image')]),
(level2conestimate, level2thresh, [('spm_mat_file',
'spm_mat_file'),
('spmT_images',
'stat_image'),
]),
(level2conestimate, datasink, [('spm_mat_file',
'2ndLevel.@spm_mat'),
('spmT_images',
'2ndLevel.@T'),
('con_images',
'2ndLevel.@con')]),
(level2thresh, datasink, [('thresholded_map',
'2ndLevel.@threshold')]),
])
And we can run it!
In [ ]:
l2analysis.run('MultiProc', plugin_args={'n_procs': 4})
Now we create a lot of outputs, but how do they look like? And also, what was the influence of different smoothing kernels and normalization?
Keep in mind, that the group analysis was only done on N=5 subjects, and that we chose a voxel-wise threshold of p<0.005. Nonetheless, we corrected for multiple comparisons with a cluster-wise FDR threshold of p<0.05.
So let's first look at the contrast congruent > incongruent:
In [ ]:
from nilearn.plotting import plot_glass_brain
plot_glass_brain(
'/output/datasink/2ndLevel/spm_con_0003_fwhm4/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='spm_fwhm4')
plot_glass_brain(
'/output/datasink/2ndLevel/spm_con_0003_fwhm8/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='spm_fwhm8')
plot_glass_brain(
'/output/datasink/2ndLevel/ants_con_0003_fwhm4/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='ants_fwhm4')
plot_glass_brain(
'/output/datasink/2ndLevel/ants_con_0003_fwhm8/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='ants_fwhm8')
Out[ ]:
The results are more or less what you would expect: The peaks are more or less at the same places for the two normalization approaches and a wider smoothing has the effect of bigger clusters, while losing the sensitivity for smaller clusters.
But if we look at the first contrast congruent, we see a different picture. In this case, the normalization with SPM seems to be more sensitive to the detection of significant voxels. Now the question is open if this is because of increased sensitivity or if this caused by an inherent normalization flaw in SPM or ANTs.
In [ ]:
from nilearn.plotting import plot_glass_brain
plot_glass_brain(
'/output/datasink/2ndLevel/spm_con_0001_fwhm4/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='spm_fwhm4')
plot_glass_brain(
'/output/datasink/2ndLevel/spm_con_0001_fwhm8/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='spm_fwhm8')
plot_glass_brain(
'/output/datasink/2ndLevel/ants_con_0001_fwhm4/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='ants_fwhm4')
plot_glass_brain(
'/output/datasink/2ndLevel/ants_con_0001_fwhm8/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='ants_fwhm8')
Out[ ]:
Last but not least, let's look at the contrast incongruent.
In [ ]:
from nilearn.plotting import plot_glass_brain
plot_glass_brain(
'/output/datasink/2ndLevel/spm_con_0002_fwhm4/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='spm_fwhm4')
plot_glass_brain(
'/output/datasink/2ndLevel/spm_con_0002_fwhm8/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='spm_fwhm8')
plot_glass_brain(
'/output/datasink/2ndLevel/ants_con_0002_fwhm4/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='ants_fwhm4')
plot_glass_brain(
'/output/datasink/2ndLevel/ants_con_0002_fwhm8/spmT_0001_thr.nii',
threshold=0, display_mode='lyrz', black_bg=True, vmax=20, title='ants_fwhm8')
Out[ ]:
Where before a smoothing of FWHM=4mm seemed to be more sensitive, here only smoothing of FWHM=8mm seems to pick up some results. Let's get some better understanding about where this significant cluster is located.
In [ ]:
from nilearn.plotting import plot_stat_map
anatimg = '/usr/share/fsl/data/standard/MNI152_T1_1mm.nii.gz'
plot_stat_map(
'/output/datasink/2ndLevel/spm_con_0002_fwhm8/spmT_0001_thr.nii'
cut_coords=(4, 6, 8, 10, 12),
bg_img=anatimg, threshold=0, display_mode='z', vmax=12)
plot_stat_map(
'/output/datasink/2ndLevel/spm_con_0002_fwhm8/spmT_0001_thr.nii'
cut_coords=(14, 16, 18, 20, 22),
bg_img=anatimg, threshold=0, display_mode='z', vmax=12)
plot_stat_map(
'/output/datasink/2ndLevel/spm_con_0002_fwhm8/spmT_0001_thr.nii',
cut_coords=(24, 26, 28, 30, 32),
bg_img=anatimg, threshold=0, display_mode='z', vmax=12)
plot_stat_map(
'/output/datasink/2ndLevel/spm_con_0002_fwhm8/spmT_0001_thr.nii',
cut_coords=(34, 36, 38, 40, 42),
bg_img=anatimg, threshold=0, display_mode='z', vmax=12)
Out[ ]: