<jupyter_start><jupyter_text>2nd-level Analysis    This notebook performes the 2nd-level analysis in template space by executing the following steps:1. Specify and estimate 2nd-level one-sample-t-test2. Threshold contrasts3. Plot output**Note:** This notebook requires that the 1st-level analysis pipeline was already executed and that it's output can be found in the dataset folder under `/dataset/derivatives/fmriflows/analysis_1stLevel/univariate`.   Data Structure RequirementsThe data structure to run this notebook should be according to the BIDS format:    dataset    ├── fmriflows_spec_analysis.json    └── derivatives        └── fmriflows            └── analysis_1stLevel                └── univariate                    └── sub-{sub_id}                        └── task-{task_id}                            └── tFilter-{tFilter_id}_sFilter-{sFilter_id}                                └── con_[con_id]_norm.nii.gz`fmriflows` will run a one-sample-t-test 2nd-level analysis on all contrasts individually.  Execution SpecificationsThis notebook will extract the relevant analysis specifications from the `fmriflows_spec_analysis.json` file in the dataset folder. In the current setup, they are as follows:<jupyter_code>import json
from os.path import join as opj

spec_file = opj('/data', 'fmriflows_spec_analysis.json')

with open(spec_file) as f:
    specs = json.load(f)
# Extract parameters for 1st-level analysis workflow
subject_list = specs['subject_list']
session_list = specs['session_list']
tasks = specs['tasks']
filters_spatial = specs['filters_spatial']
filters_temporal = specs['filters_temporal']
postfix = specs['analysis_postfix']
gm_mask_thr = specs['gm_mask_thr']
height_threshold = specs['height_threshold']
use_fwe_correction = specs['use_fwe_correction']
extent_threshold = specs['extent_threshold']
use_topo_fdr = specs['use_topo_fdr']
extent_fdr_p_threshold = specs['extent_fdr_p_threshold']
atlasreader_names = specs['atlasreader_names']
atlasreader_prob_thresh = specs['atlasreader_prob_thresh']
n_proc = specs['n_parallel_jobs']<jupyter_output><empty_output><jupyter_text>If you'd like to change any of those values manually, overwrite them below:<jupyter_code># List of subject identifiers
subject_list
# List of session identifiers
session_list
# List of spatial filters (smoothing) that were used during functional preprocessing
filters_spatial
# List of temporal filters that were used during functional preprocessing
filters_temporal
# Specify a particular analysis postfix
postfix
# Value to threshold gray matter probability template to create 2nd-level mask
gm_mask_thr
# Value for initial thresholding to define clusters
height_threshold
# Whether to use FWE (Bonferroni) correction for initial threshold
use_fwe_correction
# Minimum cluster size in voxels 
extent_threshold
# Whether to use FDR correction over cluster extent probabilities
use_topo_fdr
# P threshold to use to on FDR corrected cluster size probabilities
extent_fdr_p_threshold
# Name of atlases to use for creation of output tables
atlasreader_names
# Probability threshold to use for output tables
atlasreader_prob_thresh
# Number of parallel jobs to run
n_proc
# Task specific parameters
tasks<jupyter_output><empty_output><jupyter_text>Create the WorkflowTo ensure a good overview of the 1st-level analysis, the workflow was divided into an analysis and a report subworkflow. Import Modules<jupyter_code>from os.path import join as opj
from nipype import Node, MapNode, Workflow
from nipype.interfaces.utility import Function, IdentityInterface
from nipype.algorithms.misc import Gunzip
from nipype.interfaces.spm import OneSampleTTestDesign, EstimateModel, EstimateContrast, Threshold
from nipype.interfaces.io import SelectFiles, DataSink
# Specify SPM location
from nipype.interfaces.matlab import MatlabCommand
MatlabCommand.set_default_paths('/opt/spm12-r7219/spm12_mcr/spm12')<jupyter_output><empty_output><jupyter_text>Relevant Execution Variables<jupyter_code># Folder paths and names
exp_dir = '/data/derivatives'
out_dir = 'fmriflows'
work_dir = '/workingdir'<jupyter_output><empty_output><jupyter_text>Implement Nodes for the 2nd-level Analysis Workflow<jupyter_code># Create Mask for group analysis
def create_mask(mask_file, con_list, gm_mask_thr):

    from os.path import abspath
    from nilearn.image import resample_to_img, math_img, new_img_like
    from scipy.ndimage.morphology import binary_dilation

    # Resample mask image to contrast space and rescale to range of [0, 1]
    img_mask = resample_to_img(mask_file, con_list[0])
    mask = math_img('img/np.max(img) >= {}'.format(gm_mask_thr), img=img_mask).get_data()

    # Apply binary dilation to image
    mask = binary_dilation(mask, iterations=2)
    img_mask = new_img_like(img_mask, mask, img_mask.affine)

    # Save image as a NIfTI file
    out_file = abspath('group_mask.nii')
    img_mask.to_filename(out_file)
   
    return out_file
    
group_mask = Node(Function(input_names=['mask_file', 'con_list', 'gm_mask_thr'],
                          output_names=['out_file'],
                          function=create_mask),
                 name='group_mask')
group_mask.inputs.mask_file = '/templates/mni_icbm152_nlin_asym_09c/1.0mm_tpm_gm.nii.gz'
group_mask.inputs.gm_mask_thr = gm_mask_thr
# Gunzip NIfTI files for SPM
gunzip = MapNode(Gunzip(), name='gunzip', iterfield=['in_file'])
# Create 2nd-level desing
one_sample_ttest = Node(OneSampleTTestDesign(),
                        name="one_sample_ttest")
# Estimate 2nd-level model
level2_estimate = Node(EstimateModel(estimation_method={'Classical': 1}),
                       name="level2_estimate")
# Estimate 2nd-level contrasts
level2_con_est = Node(EstimateContrast(group_contrast=True),
                      name="level2_con_est")

cont01 = ['Group', 'T', ['mean'], [1]]
level2_con_est.inputs.contrasts = [cont01]
# Prepare contrast for two-tailed thresholding
def abs_img(spmT_file):

    from os.path import basename, abspath
    from nilearn.image import math_img

    img = math_img('np.abs(img)', img=spmT_file)
    out_file = abspath(basename(spmT_file))
    img.to_filename(out_file)

    return out_file
    
absolute_image = Node(Function(input_names=['spmT_file'],
                               output_names=['out_file'],
                               function=abs_img),
                      name='absolute_image')
# Threshold group contrast, voxel and cluster-wise
threshold = Node(Threshold(contrast_index=1,
                           use_topo_fdr=use_topo_fdr,
                           use_fwe_correction=use_fwe_correction,
                           extent_threshold=extent_threshold,
                           height_threshold=height_threshold,
                           height_threshold_type='p-value',
                           extent_fdr_p_threshold=extent_fdr_p_threshold,
                          ),
                 name='threshold')
# Mask original spmT contrast with thresholded output
def mask_img(spmT_file, thresh_image):

    from os.path import basename, abspath
    from nilearn.image import math_img

    img = math_img('img * (np.nan_to_num(thr)!=0)', img=spmT_file, thr=thresh_image)
    out_file = abspath(basename(spmT_file)).replace('.nii', '_thr.nii')
    img.to_filename(out_file)

    return out_file
    
apply_threshold = Node(Function(input_names=['spmT_file', 'thresh_image'],
                                output_names=['out_file'],
                                function=mask_img),
                       name='apply_threshold')
def apply_atlasreader(thresh_image, atlas_names, atlas_prob_thresh, extent_threshold):

    # Create output with atlasreader
    from atlasreader.atlasreader import create_output
    create_output(thresh_image,
                  cluster_extent=extent_threshold,
                  atlas=atlas_names,
                  voxel_thresh=0,
                  prob_thresh=atlas_prob_thresh)

    # Collect atlasreader output files
    from glob import glob
    out_files = glob(thresh_image.replace('_thr.nii', '_thr*.png'))
    out_files += glob(thresh_image.replace('_thr.nii', '_thr*.csv'))
    
    return out_files

atlasreader = Node(Function(input_names=['thresh_image', 'atlas_names',
                                         'atlas_prob_thresh', 'extent_threshold'],
                              output_names=['out_files'],
                              function=apply_atlasreader),
                     name='atlasreader')
atlasreader.inputs.atlas_names = atlasreader_names
atlasreader.inputs.atlas_prob_thresh = atlasreader_prob_thresh
atlasreader.inputs.extent_threshold = extent_threshold<jupyter_output><empty_output><jupyter_text>Create 2nd-level Analysis Workflow and connect nodes<jupyter_code># Create analysis workflow
analysis_2nd = Workflow(name='analysis_2nd')
analysis_2nd.base_dir = work_dir

# Add nodes to workflow and connect them
analysis_2nd.connect([(gunzip, one_sample_ttest, [('out_file', 'in_files')]),
                      (gunzip, group_mask, [('out_file', 'con_list')]),
                      (group_mask, one_sample_ttest, [('out_file', 'explicit_mask_file')]),
                      (one_sample_ttest, level2_estimate, [('spm_mat_file', 'spm_mat_file')]),
                      (level2_estimate, level2_con_est, [('spm_mat_file', 'spm_mat_file'),
                                                         ('beta_images', 'beta_images'),
                                                         ('residual_image', 'residual_image')]),
                      (level2_con_est, absolute_image, [('spmT_images', 'spmT_file')]),
                      (level2_con_est, threshold, [('spm_mat_file', 'spm_mat_file')]),
                      (absolute_image, threshold, [('out_file', 'stat_image')]),
                      (level2_con_est, apply_threshold, [('spmT_images', 'spmT_file')]),
                      (threshold, apply_threshold, [('thresholded_map', 'thresh_image')]),
                      (apply_threshold, atlasreader, [('out_file', 'thresh_image')]),
                     ])<jupyter_output><empty_output><jupyter_text>Specify Input & Output Stream<jupyter_code># Iterate over subject, session, task and run id
info_source = Node(IdentityInterface(fields=['session_id',
                                             'task_info',
                                             'spatial_filt',
                                             'temporal_filt']),
                   name='info_source')

# Generate a list of all possible contrasts, containing: (task_id, con_id, con_name)
task_info = [(t, i, c[0]) for t in list(tasks.keys())
                          for i, c in enumerate(tasks[t]['contrasts'])]

# Combine all lists of iterations
iter_list = [('task_info', task_info),
             ('spatial_filt', filters_spatial),
             ('temporal_filt', filters_temporal),
             ]

if session_list:
    iter_list.append(('session_id', session_list))
else:
    info_source.inputs.session_id = ''

info_source.iterables = iter_list
# Extract contrast specifications for 2nd-level analysis
def get_parameters(task_info):
    
    # Extract task information from (task_id, con_id)
    task_id = task_info[0]
    con_id = task_info[1] + 1

    return task_id, con_id

get_param = Node(Function(input_names=['task_info'],
                          output_names=['task_id', 'con_id'],
                          function=get_parameters),
                 name='get_param')
# Create path to input files
def create_file_path(subject_list, session_id, task_id, tFilter, sFilter, con_id, postfix):

    from glob import glob
    template_con = '/data/derivatives/fmriflows/analysis_1stLevel'
    if postfix:
        template_con += '_%s' % postfix
    template_con += '/univariate/sub-{0}/task-{1}/'
    if session_id:
        template_con += 'ses-%s/' % session_id
    template_con += '{2}_{3}/con_{4}_norm.nii???'
    
    tFilter_id = 'tFilter_%s.%s' % (tFilter[0], tFilter[1])
    sFilter_id = 'sFilter_%s.%s' % (sFilter[0], sFilter[1])

    from glob import glob
    con_files = []
    for sub_id in subject_list:
        new_cons = glob(template_con.format(
            sub_id, task_id, tFilter_id, sFilter_id, '%04d' % con_id))
        con_files += new_cons

    return sorted(con_files)

select_files = Node(Function(input_names=['subject_list', 'session_id', 'task_id',
                                          'tFilter', 'sFilter', 'con_id', 'postfix'],
                             output_names=['con_files'],
                             function=create_file_path),
                    name='select_files')
select_files.inputs.subject_list = subject_list
select_files.inputs.postfix = postfix
# Save relevant outputs in a datasink
datasink = Node(DataSink(base_directory=exp_dir,
                         container=out_dir),
                name='datasink')
# Apply the following naming substitutions for the datasink
if session_list:

    folder_old = ['_session_id_%s_spatial_filt_%s_task_info_%s_temporal_filt_%s/' % (
        ses, '.'.join([str(f) for f in sFilter]),
        '.'.join([str(t) for t in task]),
        '.'.join([str(t) for t in tFilter]))
                  for ses in session_list
                  for task in task_info
                  for sFilter in filters_spatial
                  for tFilter in filters_temporal]

    folder_new = ['task-%s/ses-%s/tFilter_%s_sFilter_%s/' % (
        '{}/{}'.format(task[0], task[2]),
        ses,
        '.'.join([str(t) for t in tFilter]),
        '.'.join([str(f) for f in sFilter]))
                  for task in task_info
                  for ses in session_list
                  for sFilter in filters_spatial
                  for tFilter in filters_temporal]
else:
    
    folder_old = ['_spatial_filt_%s_task_info_%s_temporal_filt_%s/' % (
        '.'.join([str(f) for f in sFilter]),
        '.'.join([str(t) for t in task]),
        '.'.join([str(t) for t in tFilter]))
                  for task in task_info
                  for sFilter in filters_spatial
                  for tFilter in filters_temporal]

    folder_new = ['task-%s/tFilter_%s_sFilter_%s/' % (
        '{}/{}'.format(task[0], task[2]),
        '.'.join([str(t) for t in tFilter]),
        '.'.join([str(f) for f in sFilter]))
                  for task in task_info
                  for sFilter in filters_spatial
                  for tFilter in filters_temporal]
    
substitutions = [z for z in zip(folder_old, folder_new)]
datasink.inputs.substitutions = substitutions<jupyter_output><empty_output><jupyter_text>Add Input & Output Stream to 2nd-Level Analysis Workflow<jupyter_code># Create anatomical preprocessing workflow
out_folder = 'analysis_2ndLevel'
if postfix:
    out_folder += '_%s' % postfix

# Add nodes to workflow and connect them
analysis_2nd.connect([(info_source, get_param, [('task_info', 'task_info')]),
                      (info_source, select_files, [('session_id', 'session_id'),
                                                   ('spatial_filt', 'sFilter'),
                                                   ('temporal_filt', 'tFilter')]),
                      (get_param, select_files, [('task_id', 'task_id')]),
                      (get_param, select_files, [('con_id', 'con_id')]),
                      
                      (select_files, gunzip, [('con_files', 'in_file')]),
                      
                      # Store analysis results in datasink
                      (level2_estimate, datasink, [('beta_images', '%s.univariate.@betas' % out_folder),
                                                   ('residual_image', '%s.univariate.@ResMS' % out_folder)]),
                      (level2_con_est, datasink, [('spm_mat_file', '%s.univariate.@spm_mat' % out_folder),
                                                  ('con_images', '%s.univariate.@con' % out_folder),
                                                  ('ess_images', '%s.univariate.@ess' % out_folder),
                                                  ('spmT_images', '%s.univariate.@spmT' % out_folder),
                                                  ('spmF_images', '%s.univariate.@spmF' % out_folder)]),
                      (apply_threshold, datasink, [('out_file', '%s.univariate.@thr_con' % out_folder)]),
                      (atlasreader, datasink, [('out_files', '%s.univariate.@atlas_files' % out_folder)]),
                      ])<jupyter_output><empty_output><jupyter_text>Visualize Workflow<jupyter_code># Create analysis_1st output graph
analysis_2nd.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph in the notebook
from IPython.display import Image
Image(filename=opj(analysis_2nd.base_dir, 'analysis_2nd', 'graph.png'))<jupyter_output><empty_output><jupyter_text>Run Workflow<jupyter_code># Run the workflow in parallel mode
res = analysis_2nd.run(plugin='MultiProc', plugin_args={'n_procs' : n_proc})
# Save workflow graph visualizations in datasink
analysis_2nd.write_graph(graph2use='flat', format='png', simple_form=True)
analysis_2nd.write_graph(graph2use='colored', format='png', simple_form=True)

from shutil import copyfile
copyfile(opj(analysis_2nd.base_dir, 'analysis_2nd', 'graph.png'),
         opj(exp_dir, out_dir,  out_folder, 'graph.png'))
copyfile(opj(analysis_2nd.base_dir, 'analysis_2nd', 'graph_detailed.png'),
         opj(exp_dir, out_dir, out_folder, 'graph_detailed.png'));<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10351741521539871, 'fraction_numerical': 0.005900549954170485, 'mean_word_length': 2.8485449735449735, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '474432', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5486, 'n_tokens_neox': 5223, 'n_words': 1168, 'path': 'notebooks/05_analysis_2nd-level.ipynb', 'repo_name': 'yarikoptic/fmriflows', 'size': '17456'}