from nilearn import plotting
%matplotlib inline

# Get the Node and Workflow object
from nipype import Node, Workflow

# Specify which SPM to use
from nipype.interfaces.matlab import MatlabCommand
MatlabCommand.set_default_paths('/opt/spm12-r7219/spm12_mcr/spm12')

"""
Explanation: Hands-on 2: How to create a fMRI analysis workflow
The purpose of this section is that you set up a complete fMRI analysis workflow yourself. So that in the end, you are able to perform the analysis from A-Z, i.e. from preprocessing to group analysis. This section will cover the analysis part, the previous section Hands-on 1: Preprocessing handles the preprocessing part.
We will use this opportunity to show you some nice additional interfaces/nodes that might not be relevant to your usual analysis. But it's always nice to know that they exist. And hopefully, this will encourage you to investigate all other interfaces that Nipype can bring to the tip of your finger.
<span style="color:red">Important</span>: You will not be able to go through this notebook if you haven't preprocessed your subjects first.
1st-level Analysis Workflow Structure
In this notebook, we will create a workflow that performs 1st-level analysis and normalizes the resulting beta weights to the MNI template. In concrete steps this means:
1. Specify 1st-level model parameters
2. Specify 1st-level contrasts
3. Estimate 1st-level contrasts
4. Normalize 1st-level contrasts

Imports
It's always best to have all relevant module imports at the beginning of your script. So let's import what we most certainly need.
End of explanation
"""


# Create the workflow here
# Hint: use 'base_dir' to specify where to store the working directory

analysis1st = Workflow(name='work_1st', base_dir='/output/')

"""
Explanation: Note: Ideally you would also put the imports of all the interfaces that you use here at the top. But as we will develop the workflow step by step, we can also import the relevant modules as we go.
Create Nodes and Workflow connections
Let's create all the nodes that we need! Make sure to specify all relevant inputs and keep in mind which ones you later on need to connect in your pipeline.
Workflow for the 1st-level analysis
We recommend to create the workflow and establish all its connections at a later place in your script. This helps to have everything nicely together. But for this hands-on example, it makes sense to establish the connections between the nodes as we go.
And for this, we first need to create a workflow:
End of explanation
"""


import pandas as pd
trialinfo = pd.read_table('/data/ds000114/task-fingerfootlips_events.tsv')
trialinfo

"""
Explanation: Specify 1st-level model parameters (stimuli onsets, duration, etc.)
The specify the 1st-level model we need the subject-specific onset times and duration of the stimuli. Luckily, as we are working with a BIDS dataset, this information is nicely stored in a tsv file:
End of explanation
"""


for group in trialinfo.groupby('trial_type'):
    print(group)
    print("")

"""
Explanation: Using pandas is probably the quickest and easiest ways to aggregate stimuli information per condition.
End of explanation
"""


import pandas as pd
from nipype.interfaces.base import Bunch

trialinfo = pd.read_table('/data/ds000114/task-fingerfootlips_events.tsv')
conditions = []
onsets = []
durations = []

for group in trialinfo.groupby('trial_type'):
    conditions.append(group[0])
    onsets.append(list(group[1].onset -10)) # subtracting 10s due to removing of 4 dummy scans
    durations.append(group[1].duration.tolist())

subject_info = [Bunch(conditions=conditions,
                      onsets=onsets,
                      durations=durations,
                      )]
subject_info

"""
Explanation: To create a GLM model, Nipype needs an list of Bunch objects per session. As we only have one session, our object needs to look as follows:
[Bunch(conditions=['Finger', 'Foot', 'Lips'],
       durations=[[15.0, 15.0, 15.0, 15.0, 15.0],
                  [15.0, 15.0, 15.0, 15.0, 15.0],
                  [15.0, 15.0, 15.0, 15.0, 15.0]],
       onsets=[[10, 100, 190, 280, 370],
               [40, 130, 220, 310, 400],
               [70, 160, 250, 340, 430]]
       )]

For more information see either the official documnetation or the nipype_tutorial example.
So, let's create this Bunch object that we then can use for the GLM model.
End of explanation
"""


from nipype.algorithms.modelgen import SpecifySPMModel

# Initiate the SpecifySPMModel node here

modelspec = Node(SpecifySPMModel(concatenate_runs=False,
                                 input_units='secs',
                                 output_units='secs',
                                 time_repetition=2.5,
                                 high_pass_filter_cutoff=128,
                                 subject_info=subject_info),
                 name="modelspec")

"""
Explanation: Good! Now we can create the node that will create the SPM model. For this we will be using SpecifySPMModel. As a reminder the TR of the acquisition is 2.5s and we want to use a high pass filter of 128.
End of explanation
"""


# Condition names
condition_names = ['Finger', 'Foot', 'Lips']

# Contrasts
cont01 = ['average',        'T', condition_names, [1/3., 1/3., 1/3.]]
cont02 = ['Finger',         'T', condition_names, [1, 0, 0]]
cont03 = ['Foot',           'T', condition_names, [0, 1, 0]]
cont04 = ['Lips',           'T', condition_names, [0, 0, 1]]
cont05 = ['Finger < others','T', condition_names, [-1, 0.5, 0.5]]
cont06 = ['Foot < others',  'T', condition_names, [0.5, -1, 0.5]]
cont07 = ['Lips > others',  'T', condition_names, [-0.5, -0.5, 1]]

cont08 = ['activation',     'F', [cont02, cont03, cont04]]
cont09 = ['differences',    'F', [cont05, cont06, cont07]]

contrast_list = [cont01, cont02, cont03, cont04, cont05, cont06, cont07, cont08, cont09]

"""
Explanation: This node will also need some additional inputs, such as the preprocessed functional images, the motion parameters etc. We will specify those once we take care of the workflow data input stream.
Specify 1st-level contrasts
To do any GLM analysis, we need to also define the contrasts that we want to investigate. If we recap, we had three different conditions in the fingerfootlips task in this dataset:

finger
foot
lips

Therefore, we could create the following contrasts (seven T-contrasts and two F-contrasts):
End of explanation
"""


from nipype.interfaces.spm import Level1Design

# Initiate the Level1Design node here

level1design = Node(Level1Design(bases={'hrf': {'derivs': [1, 0]}},
                                 timing_units='secs',
                                 interscan_interval=2.5,
                                 model_serial_correlations='AR(1)'),
                    name="level1design")

"""
Explanation: Estimate 1st-level contrasts
Before we can estimate the 1st-level contrasts, we first need to create the 1st-level design. Here you can also specify what kind of basis function you want (HRF, FIR, Fourier, etc.), if you want to use time and dispersion derivatives and how you want to model the serial correlation.
In this example, I propose that you use an HRF basis function, that we model time derivatives and that we model the serial correlation with AR(1).
End of explanation
"""


# Connect the two nodes here

analysis1st.connect([(modelspec, level1design, [('session_info',
                                                 'session_info')])])

"""
Explanation: Now that we have the Model Specification and 1st-Level Design node, we can connect them to each other:
End of explanation
"""


from nipype.interfaces.spm import EstimateModel

# Initiate the EstimateModel node here

level1estimate = Node(EstimateModel(estimation_method={'Classical': 1}),
                      name="level1estimate")

"""
Explanation: Now we need to estimate the model. I recommend that you'll use a Classical: 1 method to estimate the model.
End of explanation
"""


# Connect the two nodes here

analysis1st.connect([(level1design, level1estimate, [('spm_mat_file',
                                                      'spm_mat_file')])])

"""
Explanation: Now we can connect the 1st-Level Design node with the model estimation node.
End of explanation
"""


from nipype.interfaces.spm import EstimateContrast

# Initiate the EstimateContrast node here

level1conest = Node(EstimateContrast(contrasts=contrast_list),
                    name="level1conest")

"""
Explanation: Now that we estimate the model, we can estimate the contrasts. Don't forget to feed the list of contrast we specify above to this node.
End of explanation
"""


# Connect the two nodes here

analysis1st.connect([(level1estimate, level1conest, [('spm_mat_file',
                                                      'spm_mat_file'),
                                                     ('beta_images',
                                                      'beta_images'),
                                                     ('residual_image',
                                                      'residual_image')])])

"""
Explanation: Now we can connect the model estimation node with the contrast estimation node.
End of explanation
"""


from nipype.interfaces.spm import Normalize12

# Location of the template
template = '/opt/spm12-r7219/spm12_mcr/spm12/tpm/TPM.nii'

# Initiate the Normalize12 node here

normalize = Node(Normalize12(jobtype='estwrite',
                             tpm=template,
                             write_voxel_sizes=[4, 4, 4]
                            ),
                 name="normalize")

"""
Explanation: Normalize 1st-level contrasts
Now that the contrasts were estimated in subject space we can put them into a common reference space by normalizing them to a specific template. In this case, we will be using SPM12's Normalize routine and normalize to the SPM12 tissue probability map TPM.nii.
At this step, you can also specify the voxel resolution of the output volumes. If you don't specify it, it will normalize to a voxel resolution of 2x2x2mm. As a training exercise, set the voxel resolution to 4x4x4mm.
End of explanation
"""


# Connect the nodes here

analysis1st.connect([(level1conest, normalize, [('con_images',
                                                 'apply_to_files')])
                     ])

"""
Explanation: Now we can connect the estimated contrasts to normalization node.
End of explanation
"""


# Import the SelectFiles
from nipype import SelectFiles

# String template with {}-based strings
templates = {'anat': '/data/ds000114/sub-{subj_id}/ses-test/anat/sub-{subj_id}_ses-test_T1w.nii.gz',
             'func': '/output/datasink_handson/preproc/sub-{subj_id}_detrend.nii.gz',
             'mc_param': '/output/datasink_handson/preproc/sub-{subj_id}.par',
             'outliers': '/output/datasink_handson/preproc/art.sub-{subj_id}_outliers.txt'
            }

# Create SelectFiles node
sf = Node(SelectFiles(templates, sort_filelist=True),
          name='selectfiles')

"""
Explanation: Datainput with SelectFiles and iterables
As in the preprocessing hands-on, we will again be using SelectFiles and iterables. So, what do we need?
From the preprocessing pipeline, we need the functional images, the motion parameters and the list of outliers. Also, for the normalization, we need the subject-specific anatomy.
End of explanation
"""


# list of subject identifiers
subject_list = ['02', '03', '04', '07', '08', '09']
sf.iterables = [('subj_id', subject_list)]

"""
Explanation: Now we can specify over which subjects the workflow should iterate. As we preprocessed only subjects 1 to 5, we can only them for this analysis.
End of explanation
"""


from nipype.algorithms.misc import Gunzip

# Initiate the two Gunzip node here

gunzip_anat = Node(Gunzip(), name='gunzip_anat')
gunzip_func = Node(Gunzip(), name='gunzip_func')

"""
Explanation: Gunzip Node
SPM12 can accept NIfTI files as input, but online if they are not compressed ('unzipped'). Therefore, we need to use a Gunzip node to unzip the detrend file and another one to unzip the anatomy image, before we can feed it to the model specification node.
End of explanation
"""


# Connect SelectFiles node to the other nodes here

analysis1st.connect([(sf, gunzip_anat, [('anat', 'in_file')]),
                     (sf, gunzip_func, [('func', 'in_file')]),
                     (gunzip_anat, normalize, [('out_file', 'image_to_align')]),
                     (gunzip_func, modelspec, [('out_file', 'functional_runs')]),
                     (sf, modelspec, [('mc_param', 'realignment_parameters'),
                                      ('outliers', 'outlier_files'),
                                      ])
                    ])

"""
Explanation: And as a final step, we just need to connect this SelectFiles node to the rest of the workflow.
End of explanation
"""


from nipype.interfaces.io import DataSink

# Initiate DataSink node here

# Initiate the datasink node
output_folder = 'datasink_handson'
datasink = Node(DataSink(base_directory='/output/',
                         container=output_folder),
                name="datasink")

## Use the following substitutions for the DataSink output
substitutions = [('_subj_id_', 'sub-')]
datasink.inputs.substitutions = substitutions

"""
Explanation: Data output with DataSink
Now, before we run the workflow, let's again specify a Datasink folder to only keep those files that we want to keep.
End of explanation
"""


# Connect nodes to datasink here

analysis1st.connect([(level1conest, datasink, [('spm_mat_file', '1stLevel.@spm_mat'),
                                               ('spmT_images', '1stLevel.@T'),
                                               ('spmF_images', '1stLevel.@F'),
                                              ]),
                     (normalize, datasink, [('normalized_files', 'normalized.@files'),
                                            ('normalized_image', 'normalized.@image'),
                                           ]),
                    ])

"""
Explanation: Now the next step is to specify all the output that we want to keep in our output folder output. Probably best to keep are the:
- SPM.mat file and the spmT and spmF files from the contrast estimation node
- normalized betas and anatomy
End of explanation
"""


# Create 1st-level analysis output graph
analysis1st.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph
from IPython.display import Image
Image(filename='/output/work_1st/graph.png')

"""
Explanation: Visualize the workflow
Now that the workflow is finished, let's visualize it again.
End of explanation
"""


analysis1st.run('MultiProc', plugin_args={'n_procs': 4})

"""
Explanation: Run the Workflow
Now that everything is ready, we can run the 1st-level analysis workflow. Change n_procs to the number of jobs/cores you want to use.
End of explanation
"""


%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt

"""
Explanation: Visualize results
End of explanation
"""


from scipy.io import loadmat

# Using scipy's loadmat function we can access SPM.mat
spmmat = loadmat('/output/datasink_handson/1stLevel/sub-07/SPM.mat',
                 struct_as_record=False)

"""
Explanation: First, let's look at the 1st-level Design Matrix of subject one, to verify that everything is as it should be.
End of explanation
"""


designMatrix = spmmat['SPM'][0][0].xX[0][0].X
names = [i[0] for i in spmmat['SPM'][0][0].xX[0][0].name[0]]

"""
Explanation: The design matrix and the names of the regressors are a bit hidden in the spmmat variable, but they can be accessed as follows:
End of explanation
"""


normed_design = designMatrix / np.abs(designMatrix).max(axis=0)

"""
Explanation: Now before we can plot it, we just need to normalize the desing matrix in such a way, that each column has a maximum amplitude of 1. This is just for visualization purposes, otherwise the rotation parameters with their rather small values will not show up in the figure.
End of explanation
"""


fig, ax = plt.subplots(figsize=(8, 8))
plt.imshow(normed_design, aspect='auto', cmap='gray', interpolation='none')
ax.set_ylabel('Volume id')
ax.set_xticks(np.arange(len(names)))
ax.set_xticklabels(names, rotation=90);

"""
Explanation: And we're ready to plot the design matrix.
End of explanation
"""


import nibabel as nb
from nilearn.plotting import plot_anat
from nilearn.plotting import plot_glass_brain

# Load GM probability map of TPM.nii
img = nb.load('/opt/spm12-r7219/spm12_mcr/spm12/tpm/TPM.nii')
GM_template = nb.Nifti1Image(img.get_data()[..., 0], img.affine, img.header)

# Plot normalized subject anatomy
display = plot_anat('/output/datasink_handson/normalized/sub-07/wsub-07_ses-test_T1w.nii',
                    dim=-0.1)

# Overlay in edges GM map
display.add_edges(GM_template)

"""
Explanation: Now that we're happy with the design matrix, let's look how well the normalization worked.
End of explanation
"""


plot_glass_brain('/output/datasink_handson/normalized/sub-07/wess_0008.nii',
                 output_file="/output/datasink_handson/normalized/sub-07/f-contr_activation.png",
                 colorbar=True, display_mode='lyrz', black_bg=True, threshold=25,
                 title='subject 7 - F-contrast: Activation');

Image(filename='/output/datasink_handson/normalized/sub-07/f-contr_activation.png')

plot_glass_brain('/output/datasink_handson/normalized/sub-07/wess_0009.nii',
                 output_file="/output/datasink_handson/normalized/sub-07/f-contr_differences.png",
                 colorbar=True, display_mode='lyrz', black_bg=True, threshold=25,
                 title='subject 7 - F-contrast: Differences');

Image(filename='/output/datasink_handson/normalized/sub-07/f-contr_differences.png')

"""
Explanation: Let's look at the contrasts of one subject that we've just computed. In particular the F-contrast.
End of explanation
"""


# Get the Node and Workflow object
from nipype import Node, Workflow

# Specify which SPM to use
from nipype.interfaces.matlab import MatlabCommand
MatlabCommand.set_default_paths('/opt/spm12-r7219/spm12_mcr/spm12')

"""
Explanation: 2nd-level Analysis Workflow Structure
Last but not least, the group level analysis. This example will also directly include thresholding of the output, as well as some visualization.
Imports
To make sure that the necessary imports are done, here they are again:
End of explanation
"""


# Create the workflow here
# Hint: use 'base_dir' to specify where to store the working directory

analysis2nd = Workflow(name='work_2nd', base_dir='/output/')

"""
Explanation: Create Nodes and Workflow connections
Now we should know this part very well.
Workflow for the 2nd-level analysis
End of explanation
"""


from nipype.interfaces.spm import OneSampleTTestDesign

# Initiate the OneSampleTTestDesign node here

onesamplettestdes = Node(OneSampleTTestDesign(), name="onesampttestdes")

"""
Explanation: 2nd-Level Design
This step depends on your study design and the tests you want to perform. If you're using SPM to do the group analysis, you have the liberty to choose between a factorial design, a multiple regression design, one-sample T-Test design, a paired T-Test design or a two-sample T-Test design.
For the current example, we will be using a one sample T-Test design.
End of explanation
"""


from nipype.interfaces.spm import EstimateModel, EstimateContrast

# Initiate the EstimateModel and the EstimateContrast node here

level2estimate = Node(EstimateModel(estimation_method={'Classical': 1}),
                      name="level2estimate")

level2conestimate = Node(EstimateContrast(group_contrast=True),
                         name="level2conestimate")

"""
Explanation: The next two steps are the same as for the 1st-level design, i.e. estimation of the model followed by estimation of the contrasts.
End of explanation
"""


cont01 = ['Group', 'T', ['mean'], [1]]
level2conestimate.inputs.contrasts = [cont01]

"""
Explanation: To finish the EstimateContrast node, we also need to specify which contrast should be computed. For a 2nd-level one sample t-test design, this is rather straightforward:
End of explanation
"""


# Connect OneSampleTTestDesign, EstimateModel and EstimateContrast here

analysis2nd.connect([(onesamplettestdes, level2estimate, [('spm_mat_file',
                                                           'spm_mat_file')]),
                     (level2estimate, level2conestimate, [('spm_mat_file',
                                                           'spm_mat_file'),
                                                          ('beta_images',
                                                           'beta_images'),
                                                          ('residual_image',
                                                           'residual_image')])
                    ])

"""
Explanation: Now, let's connect those three design nodes to each other.
End of explanation
"""


from nipype.interfaces.spm import Threshold

level2thresh = Node(Threshold(contrast_index=1,
                              use_topo_fdr=True,
                              use_fwe_correction=False,
                              extent_threshold=0,
                              height_threshold=0.01,
                              height_threshold_type='p-value',
                              extent_fdr_p_threshold=0.05),
                    name="level2thresh")

# Connect the Threshold node to the EstimateContrast node here

analysis2nd.connect([(level2conestimate, level2thresh, [('spm_mat_file',
                                                         'spm_mat_file'),
                                                        ('spmT_images',
                                                         'stat_image'),
                                                       ])
                    ])

"""
Explanation: Thresholding of output contrast
And to close, we will use SPM Threshold. With this routine, we can set a specific voxel threshold (i.e. p<0.001) and apply an FDR cluster threshold (i.e. p<0.05).
As we only have 5 subjects, I recommend to set the voxel threshold to 0.01 and to leave the cluster threshold at 0.05.
End of explanation
"""


%%bash
TEMPLATE='/opt/spm12-r7219/spm12_mcr/spm12/tpm/TPM.nii'

# Extract the first volume with `fslroi`
fslroi $TEMPLATE GM_PM.nii.gz 0 1

# Threshold the probability mask at 10%
fslmaths GM_PM.nii -thr 0.10 -bin /output/datasink_handson/GM_mask.nii.gz

# Unzip the mask and delete the GM_PM.nii file
gunzip /output/datasink_handson/GM_mask.nii.gz
rm GM_PM.nii.gz

"""
Explanation: Gray Matter Mask
We could run our 2nd-level workflow as it is. All the major nodes are there. But I nonetheless suggest that we use a gray matter mask to restrict the analysis to only gray matter voxels.
In the 1st-level analysis, we normalized to SPM12's TPM.nii tissue probability atlas. Therefore, we could just take the gray matter probability map of this TPM.nii image (the first volume) and threshold it at a certain probability value to get a binary mask. This can of course also all be done in Nipype, but sometimes the direct bash code is quicker:
End of explanation
"""


from nilearn.plotting import plot_anat
%matplotlib inline
plot_anat('/output/datasink_handson/GM_mask.nii', dim=-1)

"""
Explanation: Let's take a look at this mask:
End of explanation
"""


onesamplettestdes.inputs.explicit_mask_file = '/output/datasink_handson/GM_mask.nii'

"""
Explanation: Now we just need to specify this binary mask as an explicit_mask_file for the one sample T-test node.
End of explanation
"""


# Import the SelectFiles
from nipype import SelectFiles

# String template with {}-based strings
templates = {'cons': '/output/datasink_handson/normalized/sub-*/w*_{cont_id}.nii'}

# Create SelectFiles node
sf = Node(SelectFiles(templates, sort_filelist=True),
          name='selectfiles')

"""
Explanation: Datainput with SelectFiles and iterables
We will again be using SelectFiles and iterables.
So, what do we need? Actually, just the 1st-level contrasts of all subjects, separated by contrast number.
End of explanation
"""


# list of contrast identifiers
contrast_id_list = ['0001', '0002', '0003', '0004', '0005',
                    '0006', '0007', '0008', '0009']
sf.iterables = [('cont_id', contrast_id_list)]

"""
Explanation: We are using * to tell SelectFiles that it can grab all available subjects and any contrast, with a specific contrast id, independnet if it's an t-contrast (con) or an F-contrast (ess) contrast.
So, let's specify over which contrast the workflow should iterate.
End of explanation
"""


analysis2nd.connect([(sf, onesamplettestdes, [('cons', 'in_files')])])

"""
Explanation: Now we need to connect the SelectFiles to the OneSampleTTestDesign node.
End of explanation
"""


from nipype.interfaces.io import DataSink

# Initiate DataSink node here

# Initiate the datasink node
output_folder = 'datasink_handson'
datasink = Node(DataSink(base_directory='/output/',
                         container=output_folder),
                name="datasink")

## Use the following substitutions for the DataSink output
substitutions = [('_cont_id_', 'con_')]
datasink.inputs.substitutions = substitutions

"""
Explanation: Data output with DataSink
Now, before we run the workflow, let's again specify a Datasink folder to only keep those files that we want to keep.
End of explanation
"""


# Connect nodes to datasink here

analysis2nd.connect([(level2conestimate, datasink, [('spm_mat_file',
                                                     '2ndLevel.@spm_mat'),
                                                    ('spmT_images',
                                                     '2ndLevel.@T'),
                                                    ('con_images',
                                                     '2ndLevel.@con')]),
                    (level2thresh, datasink, [('thresholded_map',
                                               '2ndLevel.@threshold')])
                     ])

"""
Explanation: Now the next step is to specify all the output that we want to keep in our output folder output. Probably best to keep are the:
- the SPM.mat file and the spmT images from the EstimateContrast node
- the thresholded spmT images from the Threshold node
End of explanation
"""


# Create 1st-level analysis output graph
analysis2nd.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph
from IPython.display import Image
Image(filename='/output/work_2nd/graph.png')

"""
Explanation: Visualize the workflow
And we're good to go. Let's first take a look at the workflow.
End of explanation
"""


analysis2nd.run('MultiProc', plugin_args={'n_procs': 4})

"""
Explanation: Run the Workflow
Now that everything is ready, we can run the 2nd-level analysis workflow. Change n_procs to the number of jobs/cores you want to use.
End of explanation
"""


from nilearn.plotting import plot_glass_brain
%matplotlib inline
out_path = '/output/datasink_handson/2ndLevel/'

plot_glass_brain(out_path + 'con_0001/spmT_0001_thr.nii', display_mode='lyrz',
                 black_bg=True, colorbar=True, title='average (FDR corrected)');

plot_glass_brain(out_path + 'con_0002/spmT_0001_thr.nii', display_mode='lyrz',
                 black_bg=True, colorbar=True, title='Finger (FDR corrected)');

plot_glass_brain(out_path + 'con_0003/spmT_0001_thr.nii', display_mode='lyrz',
                 black_bg=True, colorbar=True, title='Foot (FDR corrected)');

plot_glass_brain(out_path + 'con_0004/spmT_0001_thr.nii', display_mode='lyrz',
                 black_bg=True, colorbar=True, title='Lips (FDR corrected)');

plot_glass_brain(out_path + 'con_0005/spmT_0001_thr.nii', display_mode='lyrz',
                 black_bg=True, colorbar=True, title='Finger < others (FDR corrected)');

plot_glass_brain(out_path + 'con_0006/spmT_0001_thr.nii', display_mode='lyrz',
                 black_bg=True, colorbar=True, title='Foot < others (FDR corrected)');

plot_glass_brain(out_path + 'con_0007/spmT_0001_thr.nii', display_mode='lyrz',
                 black_bg=True, colorbar=True, title='Lips > others (FDR corrected)');

"""
Explanation: Visualize results
Let's take a look at the results. Keep in mind that we only have N=6 subjects and that we set the voxel threshold to a very liberal p&lt;0.01. Interpretation of the results should, therefore, be taken with a lot of caution.
End of explanation
"""


