Fixed-effects workflow for fmriprep outputs?

,
1

Hi all,

I’m writing a fixed-effects workflow that uses outputs from fmriprep. I started working from this example by @ChrisGorgolewski and re-wrote it as a Nipype workflow to make it easier to iterate over subjects and runs. As far as I could tell, the original example models a single run per session per subject; I’m having trouble adapting this code to model multiple runs per subject.

I’ve put my work so far in this repository. I’ve gotten as far as running model estimation for individual runs with FILM, but I’m not sure how to feed the outputs of FILM (filmgls) to a fixed-effects workflow (fixedfx). I built the first-level workflow so that it iterates over different runs from the start, so I think the next step is to merge the copes, varcopes, etc. across all runs of a single subject before feeding them into the fixed effects workflow.

I tried using a JoinNode (mergemodel) between filmgls and fixedfx to do this. Here’s the relevant bit from block #44:

             (mask, mergesource, [('out_file', 'mask')]),
             (filmgls, mergesource, [('copes', 'copes'),
                                    ('varcopes', 'varcopes'),
                                    ('dof_file', 'dof_file')]),
             (mergesource, mergemodel, [('mask', 'mask'),
                                        ('copes', 'copes'),
                                        ('varcopes', 'varcopes'),
                                        ('dof_file', 'dof_file')]),
             (mergemodel, fixedfx, [(('mask', pickfirst), 'flameo.mask_file'),
                                    (('copes', sort_copes), 'inputspec.copes'),
                                    ('dof_file', 'inputspec.dof_files'),
                                    (('varcopes', sort_copes), 'inputspec.varcopes'),
                                    (('copes', num_copes), 'l2model.num_copes')]),

However, I got this error:

	 [Node] Setting-up "l1.fixedfx.varcopemerge" in     "/scratch/users/nvelez/swist_cache/l1_model/task-tomloc_model-localizer_sub-04/l1/fixedfx/_run_2/varcopemerge".
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-45-5b9f69c9be41> in <module>()
----> 1 l1_workflow.run()

/src/nipype/nipype/pipeline/engine/workflows.py in run(self, plugin, plugin_args, updatehash)
    593         if str2bool(self.config['execution']['create_report']):
    594             self._write_report_info(self.base_dir, self.name, execgraph)
--> 595         runner.run(execgraph, updatehash=updatehash, config=self.config)
    596         datestr = datetime.utcnow().strftime('%Y%m%dT%H%M%S')
    597         if str2bool(self.config['execution']['write_provenance']):

/src/nipype/nipype/pipeline/plugins/linear.py in run(self, graph, config, updatehash)
     42                 if self._status_callback:
     43                     self._status_callback(node, 'start')
---> 44                 node.run(updatehash=updatehash)
     45                 if self._status_callback:
     46                     self._status_callback(node, 'end')

/src/nipype/nipype/pipeline/engine/nodes.py in run(self, updatehash)
    415         # Check hash, check whether run should be enforced
    416         logger.info('[Node] Setting-up "%s" in "%s".', self.fullname, outdir)
--> 417         cached, updated = self.is_cached()
    418 
    419         # If the node is cached, check on pklz files and finish

/src/nipype/nipype/pipeline/engine/nodes.py in is_cached(self, rm_outdated)
    301 
    302         # Update hash
--> 303         hashed_inputs, hashvalue = self._get_hashval()
    304 
    305         # The output folder does not exist: not cached

/src/nipype/nipype/pipeline/engine/nodes.py in _get_hashval(self)
   1052             return self._hashed_inputs, self._hashvalue
   1053 
-> 1054         self._check_iterfield()
   1055         hashinputs = deepcopy(self._interface.inputs)
   1056         for name in self.iterfield:

/src/nipype/nipype/pipeline/engine/nodes.py in _check_iterfield(self)
   1218             if not isdefined(getattr(self.inputs, iterfield)):
   1219                 raise ValueError(("Input %s was not set but it is listed "
-> 1220                                   "in iterfields.") % iterfield)
   1221         if len(self.iterfield) > 1:
   1222             first_len = len(

ValueError: Input in_files was not set but it is listed in iterfields.

Any suggestions on how to debug this would be much appreciated! For reference, I’ve included a diagram of the workflow at the end of the notebook.

Thanks,
Natalia

2

By looking at your code my best guess is that the function sort_copes returns an empty list. One way to debug this would be to add an assertion checking that length of the list is not zero.

1 Like
3

Probably not relevant for you anymore, but I ran into a similar issue using JoinNode (and hence found this old thread). Part of the issue is in the joinsource of the JoinNode. It should be joinsource='inputnode' in this case, not joinsource='mergesource' (the syntax here are not accurate in a couple ways). I think your script will run with these changes (code in context)

def sort_filmgls_output(copes_grouped_by_run, varcopes_grouped_by_run):
    def reshape_lists(files_grouped_by_run):
        import numpy as np
        if not isinstance(files_grouped_by_run, list):
            files = [files_grouped_by_run]
        else:
            files = files_grouped_by_run
            
        if all(len(x) == len(files[0]) for x in files): n_files = len(files[0])

        all_files = np.array(files).flatten()
        files_grouped_by_contrast = all_files.reshape(int(len(all_files) / n_files), n_files).T.tolist()
        return files_grouped_by_contrast
    
    copes_grouped_by_contrast = reshape_lists(copes_grouped_by_run)
    varcopes_grouped_by_contrast = reshape_lists(varcopes_grouped_by_run)

    return copes_grouped_by_contrast, varcopes_grouped_by_contrast

pass_run_data = pe.Node(niu.IdentityInterface(fields = ['mask', 'dof_file', 'copes', 'varcopes']), 'pass_run_data')

join_run_data = pe.JoinNode(niu.IdentityInterface(fields=['masks', 'dof_files', 'copes', 'varcopes']),  joinsource='inputnode',  joinfield=['masks', 'dof_files', 'copes', 'varcopes'], name='join_run_data')

group_by_contrast = pe.Node(niu.Function(input_names=['copes_grouped_by_run', 'varcopes_grouped_by_run'], output_names=['copes_grouped_by_contrast', 'varcopes_grouped_by_contrast'], function=sort_filmgls_output), name='group_by_contrast')

    (mask, pass_run_data, [('out_file', 'mask')]),
    (filmgls, pass_run_data, [
        ('copes', 'copes'),
        ('varcopes', 'varcopes'),
        ('dof_file', 'dof_file') ]),
    (pass_run_data, join_run_data, [
        ('mask', 'masks'),
        ('dof_file', 'dof_files'),
        ('copes', 'copes'),
        ('varcopes', 'varcopes') ]),
    (join_run_data, group_by_contrast, [
        ('copes', 'copes_grouped_by_run'),
        ('varcopes', 'varcopes_grouped_by_run')]),
    (group_by_contrast, fixed_fx, [
        ('copes_grouped_by_contrast', 'inputspec.copes'),
        ('varcopes_grouped_by_contrast', 'inputspec.varcopes')]),
    (join_run_data, fixed_fx, [
        (('masks', pickfirst), 'flameo.mask_file'),
        ('dof_files', 'inputspec.dof_files'),
        (('copes', num_copes), 'l2model.num_copes')]), ### number of runs
1 Like
4

Wow, thank you so much!! I was actually just looking at this script again this week. :slight_smile: I’ll try it and let you know how it goes.

5

HI @dae, quick side question:
You mention that the syntax here is not accurate in a couple of ways. In what way and where do you mean? Would be great if you could open either an issue or directly a PR to help us improve the tutorial.

1 Like
6

@dae: Just to follow up, your changes did work—thank you! I’m just running into a separate problem now with FLAMEO. :slight_smile:

7

Hello,

I am adapting @nataliavelez workflow to fit my own data but I’m not sure what the variable/template contrast_template is. It does not seem like a BIDS file and I’m not sure where it’s coming from. Here is the code for it, as well as where it’s used in the workflow

contrast_template = 'l1model/task-%s_model-%s.json'

datasource = pe.Node(nio.DataGrabber(infields=['subject_id',
                                               'task',
                                               'model',
                                               'run'],
                                    outfields=['struct',
                                               'func',
                                               'mask',
                                               'confounds_file',
                                               'events_file',
                                               'contrasts_file']),
                    'datasource')

datasource.inputs.base_directory = derivatives_dir
datasource.inputs.template = '*'
datasource.inputs.sort_filelist = True
datasource.inputs.field_template = dict(struct=anat_template,
                                       func=func_template,
                                       mask=mask_template,
                                       confounds_file=confounds_template,
                                       events_file=model_template,
                                       contrasts_file=contrast_template)
datasource.inputs.template_args = dict(struct=[['subject_id']],
                                      func=[['subject_id', 'task', 'run']],
                                      mask=[['subject_id', 'task', 'run']],
                                      confounds_file=[['subject_id', 'task', 'run']],
                                      events_file=[['subject_id', 'task', 'model', 'run']],
                                      contrasts_file=[['task', 'model']])

Thanks for the help.

8

@atersakyan,
If you check out the section on DataGrabber from Michael’s tutorial on basic data input, you’ll see other example usage of the datasource.inputs.field_template (about 1/3 of the way down). The explanation there echoes the documentation: field_template key, value pairs provide the iterable search paths for output fields. In this case, it seems the first level contrast estimate will be created following the convention provided by ‘contrast_template’.

Hope that points you in the right direction.

9

Hello all,

I’ve been adapting @nataliavelez script for my lab’s data, but I’ve been running into an issue with writing out the fixed effects results, where executing the code below gives me the following exeption: “Trying to connect l1.modelspc:session_info to l1.l1design:session_info but input ‘session_info’ of node ‘l1.l1design’ is already connected.” – which is fine, but the script won’t run past the exception. I only get this exception for ‘session_info’ of the l1.l1design node. When I remove the exception and force it to continue, I get a new exception telling me featmodel has no input called fsf_files. Is this something the featmodel interface should be generating? Any help is much appreciated.

(fixedfx, datasink, [
    ('outputspec.res4d', "flameo.@res4d"),
    ('outputspec.copes', "flameo.@copes"),
    ('outputspec.varcopes', "flameo.@varcopes"),
    ('output.zstats', "flameo.@zstats"),
    ('outputspec.tstats', "flameo.@tstats")
]),
10

I’ve been running into an issue with writing out the fixed effects results, where executing the code below gives me the following exeption: “Trying to connect l1.modelspc:session_info to l1.l1design:session_info but input ‘session_info’ of node ‘l1.l1design’ is already connected.”

Hmmm, I think the code that’s causing this error might be earlier in your script. Could you show us the part of your code where you’re defining the l1 workflow? In my repository, this is the chunk that starts with:

level1_workflow = pe.Workflow('l1', base_dir = work_dir)

When I remove the exception and force it to continue, I get a new exception telling me featmodel has no input called fsf_files. Is this something the featmodel interface should be generating?

What’s in featmodel.inputs? I wonder if this part of my code is just out of date. Looking through the documentation, it looks like interfaces.fsl.FEATModel has an input called “fsf_file,” not “fsf_files.”

11

Thank you for the prompt response! I figured out why I was getting the exception for already connected nodes. I simply had session_info appearing twice in that part my script, where contrasts should have been. :upside_down_face:

As for featmodel.inputs, it does look as though it calls for “fsf_file”, and not “fsf_files”. After addressing that, I no longer get that exception either. Thank you again for the help, looks like I can move forward with adapting the remainder of the script. I might reach out again if I have other questions.