diff --git a/fitlins/__about__.py b/fitlins/__about__.py index cfa465ae..16082d78 100644 --- a/fitlins/__about__.py +++ b/fitlins/__about__.py @@ -34,21 +34,20 @@ REQUIRES = [ 'nibabel>=2.0', - 'nipype>=1.1.2', + 'nipype>=1.1.6', 'seaborn>=0.7.1', 'numpy>=1.11', 'nilearn>=0.4.0', 'pandas>=0.19', 'tables>=3.2.1', - 'nistats>=0.0.1a', - 'grabbit>=0.2.3', + 'nistats>=0.0.1b0', 'pybids>=0.6.5', 'jinja2', ] LINKS_REQUIRES = [ - 'git+https://github.com/nistats/nistats.git@' - 'ce3695e8f34c6f34323766dc96a60a53b69d2729#egg=nistats', + 'git+https://github.com/bids-standard/pybids.git@' + 'dd8c9bc6478bb63e09f6c95566a11677ce12ded7#egg=pybids', ] TESTS_REQUIRES = [ diff --git a/fitlins/cli/run.py b/fitlins/cli/run.py index aa5c6959..cda961e1 100755 --- a/fitlins/cli/run.py +++ b/fitlins/cli/run.py @@ -79,9 +79,9 @@ def get_parser(): 'removed)') g_bids.add_argument('-m', '--model', action='store', default='model.json', help='location of BIDS model description (default bids_dir/model.json)') - g_bids.add_argument('-p', '--preproc-dir', action='store', default='fmriprep', - help='location of preprocessed data (if relative path, search ' - 'bids_dir/derivatives, followed by output_dir)') + g_bids.add_argument('-d', '--derivatives', action='store', nargs='+', + help='location of derivatives (including preprocessed images).' + 'If none specified, indexes all derivatives under bids_dir/derivatives.') g_bids.add_argument('--derivative-label', action='store', type=str, help='execution label to append to derivative directory name') g_bids.add_argument('--space', action='store', @@ -139,13 +139,7 @@ def run_fitlins(argv=None): if opts.model in (None, 'default') and not op.exists(model): model = 'default' - preproc_dir = default_path(opts.preproc_dir, - op.join(opts.bids_dir, 'derivatives'), - 'fmriprep') - if not op.exists(preproc_dir): - preproc_dir = default_path(opts.preproc_dir, opts.output_dir, 'fmriprep') - if not op.exists(preproc_dir): - raise RuntimeError("Preprocessed data could not be found") + derivatives = True if not opts.derivatives else opts.derivatives pipeline_name = 'fitlins' if opts.derivative_label: @@ -158,7 +152,7 @@ def run_fitlins(argv=None): work_dir = mkdtemp() if opts.work_dir is None else opts.work_dir fitlins_wf = init_fitlins_wf( - opts.bids_dir, preproc_dir, deriv_dir, opts.space, model=model, + opts.bids_dir, derivatives, deriv_dir, opts.space, model=model, participants=subject_list, base_dir=work_dir, include_pattern=opts.include, exclude_pattern=opts.exclude ) diff --git a/fitlins/interfaces/bids.py b/fitlins/interfaces/bids.py index f732d31c..0486da4a 100644 --- a/fitlins/interfaces/bids.py +++ b/fitlins/interfaces/bids.py @@ -1,5 +1,4 @@ import os -from functools import reduce from pathlib import Path from gzip import GzipFile import json @@ -7,8 +6,6 @@ import numpy as np import nibabel as nb -from collections import defaultdict - from nipype import logging from nipype.utils.filemanip import makedirs, copyfile from nipype.interfaces.base import ( @@ -18,7 +15,7 @@ ) from nipype.interfaces.io import IOBase -from ..utils import dict_intersection, snake_to_camel +from ..utils import snake_to_camel iflogger = logging.getLogger('nipype.interface') @@ -99,10 +96,11 @@ def _run_interface(self, runtime): from bids.analysis import auto_model models = self.inputs.model if not isinstance(models, list): - layout = bids.BIDSLayout(self.inputs.bids_dir) + # model is not yet standardized, so validate=False + layout = bids.BIDSLayout(self.inputs.bids_dir, validate=False) if not isdefined(models): - models = layout.get(type='model') + models = layout.get(suffix='smdl', return_type='file') if not models: raise ValueError("No models found") elif models == 'default': @@ -135,8 +133,8 @@ class LoadBIDSModelInputSpec(BaseInterfaceInputSpec): bids_dir = Directory(exists=True, mandatory=True, desc='BIDS dataset root directory') - preproc_dir = Directory(exists=True, - desc='Optional preprocessed files directory') + derivatives = traits.Either(traits.Bool, InputMultiPath(Directory(exists=True)), + desc='Derivative folders') model = traits.Dict(desc='Model specification', mandatory=True) selectors = traits.Dict(desc='Limit collected sessions', usedefault=True) include_pattern = InputMultiPath( @@ -149,8 +147,7 @@ class LoadBIDSModelInputSpec(BaseInterfaceInputSpec): class LoadBIDSModelOutputSpec(TraitedSpec): session_info = traits.List(traits.Dict()) - contrast_info = traits.List(traits.List(File())) - contrast_indices = traits.List(traits.List(traits.List(traits.Dict))) + contrast_info = traits.List(traits.List(traits.List(traits.Dict()))) entities = traits.List(traits.List(traits.Dict())) warnings = traits.List(File) @@ -161,83 +158,73 @@ class LoadBIDSModel(SimpleInterface): def _run_interface(self, runtime): import bids - bids.config.set_options(loop_preproc=True) include = self.inputs.include_pattern exclude = self.inputs.exclude_pattern + derivatives = self.inputs.derivatives if not isdefined(include): include = None if not isdefined(exclude): exclude = None + if not isdefined(derivatives): + exclude = False - paths = [(self.inputs.bids_dir, 'bids')] - if isdefined(self.inputs.preproc_dir): - paths.append((self.inputs.preproc_dir, ['bids', 'derivatives'])) - layout = bids.BIDSLayout(paths, include=include, exclude=exclude) + layout = bids.BIDSLayout(self.inputs.bids_dir, include=include, + exclude=exclude, derivatives=derivatives) selectors = self.inputs.selectors analysis = bids.Analysis(model=self.inputs.model, layout=layout) - analysis.setup(drop_na=False, **selectors) + analysis.setup(drop_na=False, desc='preproc', **selectors) self._load_level1(runtime, analysis) self._load_higher_level(runtime, analysis) - # Debug - remove, eventually - runtime.analysis = analysis - return runtime def _load_level1(self, runtime, analysis): - block = analysis.blocks[0] - block_subdir = Path(runtime.cwd) / block.level - block_subdir.mkdir(parents=True, exist_ok=True) + step = analysis.steps[0] + step_subdir = Path(runtime.cwd) / step.level + step_subdir.mkdir(parents=True, exist_ok=True) entities = [] session_info = [] - contrast_indices = [] contrast_info = [] warnings = [] - for paradigm, _, ents in block.get_design_matrix( - block.model['HRF_variables'], mode='sparse', force=True): + for sparse, dense, ents in step.get_design_matrix(): info = {} - space = analysis.layout.get_spaces(type='preproc', + space = analysis.layout.get_spaces(suffix='bold', extensions=['.nii', '.nii.gz'])[0] - preproc_files = analysis.layout.get(type='preproc', + preproc_files = analysis.layout.get(suffix='bold', extensions=['.nii', '.nii.gz'], space=space, **ents) if len(preproc_files) != 1: raise ValueError('Too many BOLD files found') - fname = preproc_files[0].filename + fname = preproc_files[0].path # Required field in seconds - TR = analysis.layout.get_metadata(fname, type='bold', + TR = analysis.layout.get_metadata(fname, suffix='bold', full_search=True)['RepetitionTime'] - dense_vars = set(block.model['variables']) - set(block.model['HRF_variables']) - - _, confounds, _ = block.get_design_matrix(dense_vars, - mode='dense', - force=True, - sampling_rate=1/TR, - **ents)[0] ent_string = '_'.join('{}-{}'.format(key, val) for key, val in ents.items()) - events_file = block_subdir / '{}_events.h5'.format(ent_string) - paradigm.to_hdf(events_file, key='events') + sparse_file = None + if sparse is not None: + sparse_file = step_subdir / '{}_sparse.h5'.format(ent_string) + sparse.to_hdf(sparse_file, key='sparse') imputed = [] - if confounds is not None: + if dense is not None: # Note that FMRIPREP includes CosineXX columns to accompany # t/aCompCor # We may want to add criteria to include HPF columns that are not # explicitly listed in the model - names = [col for col in confounds.columns - if col.startswith('NonSteadyStateOutlier') or - col in block.model['variables']] - confounds = confounds[names] + names = [col for col in dense.columns + if col.startswith('non_steady_state') or + col in step.model['x']] + dense = dense[names] # These confounds are defined pairwise with the current volume # and its predecessor, and thus may be undefined (have value @@ -246,130 +233,71 @@ def _load_level1(self, runtime, analysis): # expected NaN only. # Any other NaNs must be handled by an explicit transform in # the BIDS model. - for imputable in ('FramewiseDisplacement', - 'stdDVARS', 'non-stdDVARS', - 'vx-wisestdDVARS'): - if imputable in confounds.columns: - vals = confounds[imputable].values + for imputable in ('framewise_displacement', + 'std_dvars', 'dvars'): + if imputable in dense.columns: + vals = dense[imputable].values if not np.isnan(vals[0]): continue # Impute the mean non-zero, non-NaN value - confounds[imputable][0] = np.nanmean(vals[vals != 0]) + dense[imputable][0] = np.nanmean(vals[vals != 0]) imputed.append(imputable) - if np.isnan(confounds.values).any(): - iflogger.warning('Unexpected NaNs found in confounds; ' + if np.isnan(dense.values).any(): + iflogger.warning('Unexpected NaNs found in design matrix; ' 'regression may fail.') - confounds_file = block_subdir / '{}_confounds.h5'.format(ent_string) - confounds.to_hdf(confounds_file, key='confounds') + dense_file = step_subdir / '{}_dense.h5'.format(ent_string) + dense.to_hdf(dense_file, key='dense') else: - confounds_file = None + dense_file = None - info['events'] = str(events_file) - info['confounds'] = str(confounds_file) + info['sparse'] = str(sparse_file) if sparse_file else None + info['dense'] = str(dense_file) if dense_file else None info['repetition_time'] = TR - # Transpose so each contrast gets a row of data instead of column - contrasts, index, _ = block.get_contrasts(**ents)[0] - - contrast_type_map = defaultdict(lambda: 'T') - contrast_type_map.update({contrast['name']: contrast['type'] - for contrast in block.contrasts}) - contrast_type_list = [contrast_type_map[contrast] - for contrast in contrasts.columns] - - contrasts = contrasts.T - # Add test indicator column - contrasts['type'] = contrast_type_list + contrasts = [dict(c._asdict()) for c in step.get_contrasts(**ents)[0]] + for con in contrasts: + con['weights'] = con['weights'].to_dict('records') - contrasts_file = block_subdir / '{}_contrasts.h5'.format(ent_string) - contrasts_file.parent.mkdir(parents=True, exist_ok=True) - contrasts.to_hdf(contrasts_file, key='contrasts') - - warning_file = block_subdir / '{}_warning.html'.format(ent_string) + warning_file = step_subdir / '{}_warning.html'.format(ent_string) with warning_file.open('w') as fobj: if imputed: fobj.write(IMPUTATION_SNIPPET.format(', '.join(imputed))) entities.append(ents) session_info.append(info) - contrast_indices.append(index.to_dict('records')) - contrast_info.append(str(contrasts_file)) + contrast_info.append(contrasts) warnings.append(str(warning_file)) self._results['session_info'] = session_info self._results['warnings'] = warnings self._results.setdefault('entities', []).append(entities) - self._results.setdefault('contrast_indices', []).append(contrast_indices) self._results.setdefault('contrast_info', []).append(contrast_info) def _load_higher_level(self, runtime, analysis): - cwd = Path(runtime.cwd) - for block in analysis.blocks[1:]: - block_subdir = cwd / block.level - block_subdir.mkdir(parents=True, exist_ok=True) - - entities = [] - contrast_indices = [] + for block in analysis.steps[1:]: contrast_info = [] - for contrasts, index, ents in block.get_contrasts(): - if contrasts.empty: + for contrasts in block.get_contrasts(): + if all([c.weights.empty for c in contrasts]): continue - # The contrast index is the name of the input contrasts, - # which will very frequently be non-unique - # Hence, add the contrast to the index (table of entities) - # and switch to a matching numeric index - index['contrast'] = contrasts.index - contrasts.index = index.index - - contrast_type_map = defaultdict(lambda: 'T') - contrast_type_map.update({contrast['name']: contrast['type'] - for contrast in block.contrasts}) - contrast_type_list = [contrast_type_map[contrast] - for contrast in contrasts.columns] - - indices = index.to_dict('records') - - # Entities for a given contrast matrix include the intersection of - # entities of inputs, e.g., if this level is within-subject, the - # subject should persist - out_ents = reduce(dict_intersection, indices) - # Explicit entities take precedence over derived - out_ents.update(ents) - # Input-level contrasts will be overridden by the current level - out_ents.pop('contrast', None) - - ent_string = '_'.join('{}-{}'.format(key, val) - for key, val in out_ents.items()) - - # Transpose so each contrast gets a row of data instead of column - contrasts = contrasts.T - # Add test indicator column - contrasts['type'] = contrast_type_list - - contrasts_file = block_subdir / '{}_contrasts.h5'.format(ent_string) - contrasts_file.parent.mkdir(parents=True, exist_ok=True) - contrasts.to_hdf(contrasts_file, key='contrasts') - - entities.append(out_ents) - contrast_indices.append(indices) - contrast_info.append(str(contrasts_file)) - - self._results['entities'].append(entities) + contrasts = [dict(c._asdict()) for c in contrasts] + for contrast in contrasts: + contrast['weights'] = contrast['weights'].to_dict('records') + contrast_info.append(contrasts) + self._results['contrast_info'].append(contrast_info) - self._results['contrast_indices'].append(contrast_indices) class BIDSSelectInputSpec(BaseInterfaceInputSpec): bids_dir = Directory(exists=True, mandatory=True, desc='BIDS dataset root directories') - preproc_dir = Directory(exists=True, - desc='Optional preprocessed files directory') + derivatives = traits.Either(True, InputMultiPath(Directory(exists=True)), + desc='Derivative folders') entities = InputMultiPath(traits.Dict(), mandatory=True) selectors = traits.Dict(desc='Additional selectors to be applied', usedefault=True) @@ -387,10 +315,9 @@ class BIDSSelect(SimpleInterface): def _run_interface(self, runtime): import bids - paths = [(self.inputs.bids_dir, 'bids')] - if isdefined(self.inputs.preproc_dir): - paths.append((self.inputs.preproc_dir, ['bids', 'derivatives'])) - layout = bids.BIDSLayout(paths) + + derivatives = self.inputs.derivatives + layout = bids.BIDSLayout(self.inputs.bids_dir, derivatives=derivatives) bold_files = [] mask_files = [] @@ -410,19 +337,20 @@ def _run_interface(self, runtime): "".format(self.inputs.bids_dir, selectors, "\n\t".join( '{} ({})'.format( - f.filename, - layout.files[f.filename].entities) + f.path, + layout.files[f.path].entities) for f in bold_file))) # Select exactly matching mask file (may be over-cautious) - bold_ents = layout.parse_file_entities( - bold_file[0].filename) - bold_ents['type'] = 'brainmask' + bold_ents = layout.parse_file_entities(bold_file[0].path) + bold_ents['suffix'] = 'mask' + bold_ents['desc'] = 'brain' mask_file = layout.get(extensions=['.nii', '.nii.gz'], **bold_ents) - bold_ents.pop('type') + bold_ents.pop('suffix') + bold_ents.pop('desc') - bold_files.append(bold_file[0].filename) - mask_files.append(mask_file[0].filename if mask_file else None) + bold_files.append(bold_file[0].path) + mask_files.append(mask_file[0].path if mask_file else None) entities.append(bold_ents) self._results['bold_files'] = bold_files @@ -490,7 +418,7 @@ def _list_outputs(self): os.makedirs(base_dir, exist_ok=True) - layout = bids.BIDSLayout(base_dir) + layout = bids.BIDSLayout(base_dir, validate=False) path_patterns = self.inputs.path_patterns if not isdefined(path_patterns): path_patterns = None diff --git a/fitlins/interfaces/nistats.py b/fitlins/interfaces/nistats.py index 883a0857..1065f6e8 100644 --- a/fitlins/interfaces/nistats.py +++ b/fitlins/interfaces/nistats.py @@ -1,5 +1,4 @@ import os -from functools import reduce import numpy as np import pandas as pd import nibabel as nb @@ -9,80 +8,45 @@ from nipype.interfaces.base import ( LibraryBaseInterface, SimpleInterface, BaseInterfaceInputSpec, TraitedSpec, - OutputMultiObject, File, traits, isdefined + File, traits, isdefined ) -from ..utils import dict_intersection - class NistatsBaseInterface(LibraryBaseInterface): _pkg = 'nistats' -def build_contrast_matrix(contrast_spec, design_matrix, - identity=None): - """Construct contrast matrix and return contrast type - - Parameters - ---------- - contrast_spec : DataFrame - Weight matrix with contrasts as rows and regressors as columns - May have 'type' column indicating T/F test - design_matrix : DataFrame - GLM design matrix with regressors as columns and TR time as rows - identity : list of strings - Names of explanatory variables to ensure "identity" contrasts are - provided. - - Returns - ------- - contrast_matrix : DataFrame - Weight matrix with contrasts as columns and regressors as rows. - Regressors match columns (including order) of design matrix. - Identity contrasts are included. - contrast_types : Series - Series of 'T'/'F' indicating the type of test for each column in - the returned contrast matrix. - Identity contrasts use T-tests. +def prepare_contrasts(contrasts, all_regressors): + """ Make mutable copy of contrast list, and + generate contrast design_matrix from dictionary weight mapping """ - # The basic spec is just a transposed matrix with an optional 'type' - # column - # We'll re-transpose and expand this matrix - init_matrix = contrast_spec.drop('type', axis='columns').T - init_types = contrast_spec['type'] if 'type' in contrast_spec \ - else pd.Series() - - if identity is None: - identity = [] - all_cols = init_matrix.columns.tolist() - all_cols.extend(set(identity) - set(all_cols)) - - contrast_matrix = pd.DataFrame(index=design_matrix.columns, - columns=all_cols, data=0) - contrast_matrix.loc[tuple(init_matrix.axes)] = init_matrix - - contrast_types = pd.Series(index=all_cols, data='T') - contrast_types[init_types.index] = init_types + if not isdefined(contrasts): + return [] - if identity: - contrast_matrix.loc[identity, identity] = np.eye(len(identity)) - contrast_types[identity] = 'T' + out_contrasts = [] + for contrast in contrasts: + # Fill in zeros + weights = np.array([ + [row[col] if col in row else 0 for col in all_regressors] + for row in contrast['weights'] + ]) + out_contrasts.append( + (contrast['name'], weights, contrast['type'])) - return contrast_matrix, contrast_types + return out_contrasts class FirstLevelModelInputSpec(BaseInterfaceInputSpec): bold_file = File(exists=True, mandatory=True) mask_file = File(exists=True) session_info = traits.Dict() - contrast_info = File(exists=True) + contrast_info = traits.List(traits.Dict) class FirstLevelModelOutputSpec(TraitedSpec): contrast_maps = traits.List(File) contrast_metadata = traits.List(traits.Dict) design_matrix = File() - contrast_matrix = File() class FirstLevelModel(NistatsBaseInterface, SimpleInterface): @@ -91,50 +55,38 @@ class FirstLevelModel(NistatsBaseInterface, SimpleInterface): def _run_interface(self, runtime): info = self.inputs.session_info - img = nb.load(self.inputs.bold_file) vols = img.shape[3] - events = pd.read_hdf(info['events'], key='events') - - if info['confounds'] is not None and info['confounds'] != 'None': - confounds = pd.read_hdf(info['confounds'], key='confounds') - confound_names = confounds.columns.tolist() - drift_model = None if 'Cosine00' in confound_names else 'cosine' + if info['sparse'] not in (None, 'None'): + sparse = pd.read_hdf(info['sparse'], key='sparse').rename( + columns={'condition': 'trial_type', + 'amplitude': 'modulation'}) + sparse = sparse.dropna(subset=['modulation']) # Drop NAs else: - confounds = None - confound_names = None - drift_model = 'cosine' + sparse = None - if isdefined(self.inputs.contrast_info): - contrast_spec = pd.read_hdf(self.inputs.contrast_info, - key='contrasts') + if info['dense'] not in (None, 'None'): + dense = pd.read_hdf(info['dense'], key='dense') + column_names = dense.columns.tolist() + drift_model = None if 'cosine_00' in column_names else 'cosine' else: - contrast_spec = pd.DataFrame() + dense = None + column_names = None + drift_model = 'cosine' - mat = dm.make_design_matrix( + mat = dm.make_first_level_design_matrix( frame_times=np.arange(vols) * info['repetition_time'], - paradigm=events.rename(columns={'condition': 'trial_type', - 'amplitude': 'modulation'}), - add_regs=confounds, - add_reg_names=confound_names, + events=sparse, + add_regs=dense, + add_reg_names=column_names, drift_model=drift_model, - ) - - # Assume that explanatory variables == HRF-convolved variables - exp_vars = events['condition'].unique().tolist() - - contrast_matrix, contrast_types = build_contrast_matrix(contrast_spec, - mat, exp_vars) + ) mat.to_csv('design.tsv', sep='\t') self._results['design_matrix'] = os.path.join(runtime.cwd, 'design.tsv') - contrast_matrix.to_csv('contrasts.tsv', sep='\t') - self._results['contrast_matrix'] = os.path.join( - runtime.cwd, 'contrasts.tsv') - mask_file = self.inputs.mask_file if not isdefined(mask_file): mask_file = None @@ -143,17 +95,22 @@ def _run_interface(self, runtime): contrast_maps = [] contrast_metadata = [] - stat_fmt = os.path.join(runtime.cwd, '{}.nii.gz').format - for contrast, ctype in zip(contrast_matrix, contrast_types): - es = flm.compute_contrast(contrast_matrix[contrast].values, - {'T': 't', 'F': 'F'}[ctype], - output_type='effect_size') - es_fname = stat_fmt(contrast) + out_ents = self.inputs.contrast_info[0]['entities'] + for name, weights, type in prepare_contrasts( + self.inputs.contrast_info, mat.columns.tolist()): + es = flm.compute_contrast( + weights, type, output_type='effect_size') + es_fname = os.path.join( + runtime.cwd, '{}.nii.gz').format(name) es.to_filename(es_fname) contrast_maps.append(es_fname) - contrast_metadata.append({'contrast': contrast, - 'type': 'effect'}) + contrast_metadata.append( + {'contrast': name, + 'type': 'effect', + **out_ents} + ) + self._results['contrast_maps'] = contrast_maps self._results['contrast_metadata'] = contrast_metadata @@ -162,9 +119,8 @@ def _run_interface(self, runtime): class SecondLevelModelInputSpec(BaseInterfaceInputSpec): stat_files = traits.List(traits.List(File(exists=True)), mandatory=True) - stat_metadata = traits.List(traits.List(traits.Dict)) - contrast_info = File(exists=True) - contrast_indices = traits.List(traits.Dict) + stat_metadata = traits.List(traits.List(traits.Dict), mandatory=True) + contrast_info = traits.List(traits.Dict, mandatory=True) class SecondLevelModelOutputSpec(TraitedSpec): @@ -190,55 +146,37 @@ class SecondLevelModel(NistatsBaseInterface, SimpleInterface): def _run_interface(self, runtime): model = level2.SecondLevelModel() - files = [] - # Super inefficient... think more about this later - for idx in self.inputs.contrast_indices: - for fname, metadata in zip(_flatten(self.inputs.stat_files), - _flatten(self.inputs.stat_metadata)): - if _match(idx, metadata): - files.append(fname) - break - else: - raise ValueError - - out_ents = reduce(dict_intersection, self.inputs.contrast_indices) - in_ents = [{key: val for key, val in index.items() if key not in out_ents} - for index in self.inputs.contrast_indices] - - contrast_spec = pd.read_hdf(self.inputs.contrast_info, - key='contrasts') - - contrast_matrix = contrast_spec.drop(columns=['type']).T - contrast_types = contrast_spec['type'] - - contrast_matrix.index = ['_'.join('{}-{}'.format(key, val) - for key, val in ents.items()) - for ents in in_ents] - contrast_matrix.to_csv('contrasts.tsv', sep='\t') - self._results['contrast_matrix'] = os.path.join( - runtime.cwd, 'contrasts.tsv') - - out_ents['type'] = 'stat' - contrast_maps = [] contrast_metadata = [] - stat_fmt = os.path.join(runtime.cwd, '{}.nii.gz').format - for contrast, ctype in zip(contrast_matrix, contrast_types): - intercept = contrast_matrix[contrast] - data = np.array(files)[intercept != 0].tolist() - intercept = intercept[intercept != 0] - - model.fit(data, design_matrix=pd.DataFrame({'intercept': intercept})) - stat_type = {'T': 't', 'F': 'F'}[ctype] - stat = model.compute_contrast(second_level_stat_type=stat_type) - stat_fname = stat_fmt(contrast) + entities = self.inputs.contrast_info[0]['entities'] # Same for all + out_ents = {'type': 'stat', **entities} + + # Only keep files which match all entities for contrast + stat_metadata = _flatten(self.inputs.stat_metadata) + stat_files = _flatten(self.inputs.stat_files) + filtered_files = [] + names = [] + for m, f in zip(stat_metadata, stat_files): + if _match(entities, m): + filtered_files.append(f) + names.append(m['contrast']) + + for name, weights, type in prepare_contrasts(self.inputs.contrast_info, names): + # Need to add F-test support for intercept (more than one column) + # Currently only taking 0th column as intercept (t-test) + weights = weights[0] + input = (np.array(filtered_files)[weights != 0]).tolist() + design_matrix = pd.DataFrame({'intercept': weights[weights != 0]}) + + model.fit(input, design_matrix=design_matrix) + + stat = model.compute_contrast(second_level_stat_type=type) + stat_fname = os.path.join(runtime.cwd, '{}.nii.gz').format(name) stat.to_filename(stat_fname) contrast_maps.append(stat_fname) - metadata = out_ents.copy() - metadata['contrast'] = contrast - contrast_metadata.append(metadata) + contrast_metadata.append({'contrast': name, **out_ents}) self._results['contrast_maps'] = contrast_maps self._results['contrast_metadata'] = contrast_metadata diff --git a/fitlins/viz/reports.py b/fitlins/viz/reports.py index d686a000..a874b915 100644 --- a/fitlins/viz/reports.py +++ b/fitlins/viz/reports.py @@ -2,7 +2,7 @@ from pathlib import Path import jinja2 import pkg_resources as pkgr -import bids +from bids.layout import add_config_paths, BIDSLayout from ..utils import snake_to_camel @@ -11,6 +11,8 @@ 'model-{model}.html' ] +add_config_paths(fitlins=pkgr.resource_filename('fitlins', 'data/fitlins.json')) + def deroot(val, root): if isinstance(val, str): @@ -28,10 +30,12 @@ def deroot(val, root): def parse_directory(deriv_dir, work_dir, analysis): - fl_layout = bids.BIDSLayout( - (deriv_dir, ['bids', 'derivatives', - pkgr.resource_filename('fitlins', 'data/fitlins.json')])) - wd_layout = bids.BIDSLayout(str(Path(work_dir) / 'reportlets' / 'fitlins')) + fl_layout = BIDSLayout( + deriv_dir, + config=['bids', 'derivatives', 'fitlins']) + wd_layout = BIDSLayout( + str(Path(work_dir) / 'reportlets' / 'fitlins'), + validate=False) contrast_svgs = fl_layout.get(extensions='.svg', type='contrasts') analyses = [] @@ -73,9 +77,8 @@ def parse_directory(deriv_dir, work_dir, analysis): def write_report(level, report_dicts, run_context, deriv_dir): - fl_layout = bids.BIDSLayout( - (deriv_dir, ['bids', 'derivatives', - pkgr.resource_filename('fitlins', 'data/fitlins.json')])) + fl_layout = BIDSLayout( + deriv_dir, config=['bids', 'derivatives', 'fitlins']) env = jinja2.Environment( loader=jinja2.FileSystemLoader( diff --git a/fitlins/workflows/base.py b/fitlins/workflows/base.py index afdf90f8..577b8ef7 100644 --- a/fitlins/workflows/base.py +++ b/fitlins/workflows/base.py @@ -9,7 +9,7 @@ from ..interfaces.utils import MergeAll -def init_fitlins_wf(bids_dir, preproc_dir, out_dir, space, exclude_pattern=None, +def init_fitlins_wf(bids_dir, derivatives, out_dir, space, exclude_pattern=None, include_pattern=None, model=None, participants=None, base_dir=None, name='fitlins_wf'): wf = pe.Workflow(name=name, base_dir=base_dir) @@ -22,6 +22,9 @@ def init_fitlins_wf(bids_dir, preproc_dir, out_dir, space, exclude_pattern=None, all_models = specs.run().outputs.model_spec if not all_models: raise RuntimeError("Unable to find or construct models") + if isinstance(all_models, list): + raise RuntimeError("Currently unable to run multiple models in parallel - " + "please specify model") # # Load and run the model @@ -31,35 +34,26 @@ def init_fitlins_wf(bids_dir, preproc_dir, out_dir, space, exclude_pattern=None, loader = pe.Node( LoadBIDSModel(bids_dir=bids_dir, - preproc_dir=preproc_dir, - selectors=selectors), + derivatives=derivatives, + selectors=selectors, + model=all_models), name='loader') - if preproc_dir is not None: - loader.inputs.preproc_dir = preproc_dir if exclude_pattern is not None: loader.inputs.exclude_pattern = exclude_pattern if include_pattern is not None: loader.inputs.include_pattern = include_pattern - if isinstance(all_models, list): - loader.iterables = ('model', all_models) - else: - loader.inputs.model = all_models - # Because pybids generates the entire model in one go, we will need # various helper nodes to select the correct portions of the model select_l1_entities = pe.Node(niu.Select(index=0), name='select_l1_entities') # Select preprocessed BOLD series to analyze getter = pe.Node( - BIDSSelect(bids_dir=bids_dir, - selectors={'type': 'preproc', 'space': space}), + BIDSSelect(bids_dir=bids_dir, derivatives=derivatives, + selectors={'type': 'preproc', 'suffix': 'bold', 'space': space}), name='getter') - if preproc_dir is not None: - getter.inputs.preproc_dir = preproc_dir - select_l1_contrasts = pe.Node(niu.Select(index=0), name='select_l1_contrasts') # Run first level model @@ -82,15 +76,13 @@ def join_dict(base_dict, dict_list): collate_first_level = pe.Node(MergeAll(['contrast_maps', 'contrast_metadata']), name='collate_first_level') - select_l2_entities = pe.Node(niu.Select(index=1), name='select_l2_entities') - select_l2_indices = pe.Node(niu.Select(index=1), name='select_l2_indices') select_l2_contrasts = pe.Node(niu.Select(index=1), name='select_l2_contrasts') # Run second-level model # TODO: Iterate over all higher levels l2_model = pe.MapNode( SecondLevelModel(), - iterfield=['contrast_info', 'contrast_indices'], + iterfield=['contrast_info'], name='l2_model') collate_second_level = pe.Node(MergeAll(['contrast_maps', 'contrast_metadata']), @@ -105,34 +97,36 @@ def join_dict(base_dict, dict_list): iterfield='data', name='plot_design') - def _get_evs(info): - import pandas as pd - events = pd.read_hdf(info['events'], key='events') - return len(events['condition'].unique()) + # TODO: get_evs should be based on contrasts, EVs are any used in contrasts, + # others are confounds + # def _get_evs(info): + # import pandas as pd + # events = pd.read_hdf(info['events'], key='events') + # return len(events['condition'].unique()) - # Number of explanatory variables is used to mark off sections of the - # correlation matrix - get_evs = pe.MapNode(niu.Function(function=_get_evs), iterfield='info', name='get_evs') + # # Number of explanatory variables is used to mark off sections of the + # # correlation matrix + # get_evs = pe.MapNode(niu.Function(function=_get_evs), iterfield='info', name='get_evs') - plot_corr = pe.MapNode( - DesignCorrelationPlot(image_type='svg'), - iterfield=['data', 'explanatory_variables'], - name='plot_corr') + # plot_corr = pe.MapNode( + # DesignCorrelationPlot(image_type='svg'), + # iterfield=['data', 'explanatory_variables'], + # name='plot_corr') - plot_l1_contrast_matrix = pe.MapNode( - ContrastMatrixPlot(image_type='svg'), - iterfield='data', - name='plot_l1_contrast_matrix') + # plot_l1_contrast_matrix = pe.MapNode( + # ContrastMatrixPlot(image_type='svg'), + # iterfield='data', + # name='plot_l1_contrast_matrix') plot_l1_contrasts = pe.MapNode( GlassBrainPlot(image_type='png'), iterfield='data', name='plot_l1_contrasts') - plot_l2_contrast_matrix = pe.MapNode( - ContrastMatrixPlot(image_type='svg'), - iterfield='data', - name='plot_l2_contrast_matrix') + # plot_l2_contrast_matrix = pe.MapNode( + # ContrastMatrixPlot(image_type='svg'), + # iterfield='data', + # name='plot_l2_contrast_matrix') plot_l2_contrasts = pe.MapNode( GlassBrainPlot(image_type='png'), @@ -146,7 +140,6 @@ def _get_evs(info): reportlet_dir = Path(base_dir) / 'reportlets' / 'fitlins' reportlet_dir.mkdir(parents=True, exist_ok=True) - snippet_pattern = '[sub-{subject}/][ses-{session}/][sub-{subject}_]' \ '[ses-{session}_]task-{task}_[run-{run}_]snippet.html' ds_model_warnings = pe.MapNode( @@ -189,28 +182,28 @@ def _get_evs(info): run_without_submitting=True, name='ds_design') - ds_corr = pe.MapNode( - BIDSDataSink(base_directory=out_dir, fixed_entities={'type': 'corr'}, - path_patterns=image_pattern), - iterfield=['entities', 'in_file'], - run_without_submitting=True, - name='ds_corr') - - ds_l1_contrasts = pe.MapNode( - BIDSDataSink(base_directory=out_dir, - fixed_entities={'type': 'contrasts'}, - path_patterns=image_pattern), - iterfield=['entities', 'in_file'], - run_without_submitting=True, - name='ds_l1_contrasts') - - ds_l2_contrasts = pe.MapNode( - BIDSDataSink(base_directory=out_dir, - fixed_entities={'type': 'contrasts'}, - path_patterns=image_pattern), - iterfield=['entities', 'in_file'], - run_without_submitting=True, - name='ds_l2_contrasts') + # ds_corr = pe.MapNode( + # BIDSDataSink(base_directory=out_dir, fixed_entities={'type': 'corr'}, + # path_patterns=image_pattern), + # iterfield=['entities', 'in_file'], + # run_without_submitting=True, + # name='ds_corr') + + # ds_l1_contrasts = pe.MapNode( + # BIDSDataSink(base_directory=out_dir, + # fixed_entities={'type': 'contrasts'}, + # path_patterns=image_pattern), + # iterfield=['entities', 'in_file'], + # run_without_submitting=True, + # name='ds_l1_contrasts') + + # ds_l2_contrasts = pe.MapNode( + # BIDSDataSink(base_directory=out_dir, + # fixed_entities={'type': 'contrasts'}, + # path_patterns=image_pattern), + # iterfield=['entities', 'in_file'], + # run_without_submitting=True, + # name='ds_l2_contrasts') contrast_plot_pattern = '[sub-{subject}/][ses-{session}/]' \ '[sub-{subject}_][ses-{session}_]task-{task}[_acq-{acquisition}]' \ @@ -253,13 +246,10 @@ def _get_evs(info): (l1_model, collate_first_level, [('contrast_maps', 'contrast_maps')]), (l1_metadata, collate_first_level, [('out', 'contrast_metadata')]), - (loader, select_l2_entities, [('entities', 'inlist')]), - (loader, select_l2_indices, [('contrast_indices', 'inlist')]), (loader, select_l2_contrasts, [('contrast_info', 'inlist')]), (l1_model, l2_model, [('contrast_maps', 'stat_files')]), (l1_metadata, l2_model, [('out', 'stat_metadata')]), - (select_l2_indices, l2_model, [('out', 'contrast_indices')]), (select_l2_contrasts, l2_model, [('out', 'contrast_info')]), (l2_model, collate_second_level, [('contrast_maps', 'contrast_maps'), @@ -270,15 +260,15 @@ def _get_evs(info): # (l1_model, plot_design, [('design_matrix', 'data')]), - (loader, get_evs, [('session_info', 'info')]), - (l1_model, plot_corr, [('design_matrix', 'data')]), - (get_evs, plot_corr, [('out', 'explanatory_variables')]), + # (loader, get_evs, [('session_info', 'info')]), + # (l1_model, plot_corr, [('design_matrix', 'data')]), + # (get_evs, plot_corr, [('out', 'explanatory_variables')]), - (l1_model, plot_l1_contrast_matrix, [('contrast_matrix', 'data')]), + # (l1_model, plot_l1_contrast_matrix, [('contrast_matrix', 'data')]), (collate_first_level, plot_l1_contrasts, [('contrast_maps', 'data')]), - (l2_model, plot_l2_contrast_matrix, [('contrast_matrix', 'data')]), + # (l2_model, plot_l2_contrast_matrix, [('contrast_matrix', 'data')]), (collate_second_level, plot_l2_contrasts, [('contrast_maps', 'data')]), @@ -300,18 +290,17 @@ def _get_evs(info): (select_l1_entities, ds_design, [('out', 'entities')]), (plot_design, ds_design, [('figure', 'in_file')]), - (select_l1_entities, ds_corr, [('out', 'entities')]), - (plot_corr, ds_corr, [('figure', 'in_file')]), + # (select_l1_entities, ds_corr, [('out', 'entities')]), + # (plot_corr, ds_corr, [('figure', 'in_file')]), - (select_l1_entities, ds_l1_contrasts, [('out', 'entities')]), - (plot_l1_contrast_matrix, ds_l1_contrasts, [('figure', 'in_file')]), + # (select_l1_entities, ds_l1_contrasts, [('out', 'entities')]), + # (plot_l1_contrast_matrix, ds_l1_contrasts, [('figure', 'in_file')]), (collate_first_level, ds_l1_contrast_plots, [('contrast_metadata', 'entities')]), (plot_l1_contrasts, ds_l1_contrast_plots, [('figure', 'in_file')]), - (select_l2_entities, ds_l2_contrasts, [('out', 'entities')]), - (plot_l2_contrast_matrix, ds_l2_contrasts, [('figure', 'in_file')]), + # (plot_l2_contrast_matrix, ds_l2_contrasts, [('figure', 'in_file')]), (collate_second_level, ds_l2_contrast_plots, [('contrast_metadata', 'entities')]), (plot_l2_contrasts, ds_l2_contrast_plots, [('figure', 'in_file')]), diff --git a/requirements.txt b/requirements.txt index 330dfb98..91570e63 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,7 +3,6 @@ seaborn>=0.7.1 numpy>=1.11 nilearn>=0.4.0 pandas>=0.19 -nipype>=1.1.2 -grabbit>=0.2.3 -pybids>=0.6.5 -git+https://github.com/nistats/nistats.git@ce3695e8f34c6f34323766dc96a60a53b69d2729#egg=nistats +nipype>=1.1.6 +git+https://github.com/nistats/nistats.git@009ce3fddb3dd01e82bca3ad7d2cdbeece0138f2#egg=nistats +git+https://github.com/bids-standard/pybids.git@dd8c9bc6478bb63e09f6c95566a11677ce12ded7#egg=pybids