From 80e7adbbe233cad9a3d945563a55d786df4a3d52 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 3 Sep 2024 16:11:09 -0400 Subject: [PATCH] Update. --- env.yml | 3 - scripts/fetch_templates.py | 4 +- src/fmripost_rapidtide/cli/parser.py | 22 ++--- src/fmripost_rapidtide/config.py | 22 ++--- src/fmripost_rapidtide/data/io_spec.json | 38 ++++++--- src/fmripost_rapidtide/interfaces/nilearn.py | 16 +++- .../interfaces/rapidtide.py | 24 ++++-- src/fmripost_rapidtide/workflows/base.py | 3 + src/fmripost_rapidtide/workflows/rapidtide.py | 80 ++++--------------- 9 files changed, 99 insertions(+), 113 deletions(-) diff --git a/env.yml b/env.yml index 7924ee5..e572616 100644 --- a/env.yml +++ b/env.yml @@ -23,9 +23,6 @@ dependencies: # Utilities - graphviz=9.0 - pandoc=3.1 - # Workflow dependencies: FSL (versions pinned in 6.0.7.7) - - fsl-bet2=2111.4 - - fsl-miscmaths=2203.2 # Workflow dependencies: ANTs - ants=2.5 # Other dependencies diff --git a/scripts/fetch_templates.py b/scripts/fetch_templates.py index d622d8a..c1b5c2b 100755 --- a/scripts/fetch_templates.py +++ b/scripts/fetch_templates.py @@ -44,10 +44,10 @@ def fetch_MNI6(): """ template = 'MNI152NLin6Asym' - tf.get(template, resolution=(1, 2), desc=None, suffix='T1w') + tf.get(template, resolution=(1, 2), desc=[None, 'brain'], suffix='T1w') tf.get(template, resolution=(1, 2), desc='brain', suffix='mask') # CIFTI - tf.get(template, resolution=2, atlas='HCP', suffix='dseg') + # tf.get(template, resolution=2, atlas='HCP', suffix='dseg') # Transform from MNI152NLin2009cAsym to MNI152NLin6Asym tf.get( template, mode='image', suffix='xfm', extension='.h5', **{'from': 'MNI152NLin2009cAsym'} diff --git a/src/fmripost_rapidtide/cli/parser.py b/src/fmripost_rapidtide/cli/parser.py index 0f803a5..399d189 100644 --- a/src/fmripost_rapidtide/cli/parser.py +++ b/src/fmripost_rapidtide/cli/parser.py @@ -235,7 +235,7 @@ def _bids_filter(value, parser): ) filt_opts.add_argument( '--filterfreqs', - dest='passvec', + dest='filterfreqs', action='store', nargs=2, type=float, @@ -249,7 +249,7 @@ def _bids_filter(value, parser): ) filt_opts.add_argument( '--filterstopfreqs', - dest='stopvec', + dest='filterstopfreqs', action='store', nargs=2, type=float, @@ -266,7 +266,7 @@ def _bids_filter(value, parser): sigcalc_opts = g_rapidtide.add_argument_group('Significance calculation options') sigcalc_opts.add_argument( '--numnull', - dest='numestreps', + dest='numnull', action='store', type=int, metavar='NREPS', @@ -288,7 +288,7 @@ def _bids_filter(value, parser): ) preproc.add_argument( '--spatialfilt', - dest='gausssigma', + dest='spatialfilt', action='store', type=float, metavar='GAUSSSIGMA', @@ -302,7 +302,7 @@ def _bids_filter(value, parser): ) preproc.add_argument( '--confoundfile', - dest='confoundfilespec', + dest='confoundfile', metavar='CONFFILE', help=( 'Read additional (non-motion) confound regressors out of CONFFILE file ' @@ -314,7 +314,7 @@ def _bids_filter(value, parser): ) preproc.add_argument( '--confoundpowers', - dest='confound_power', + dest='confoundpowers', metavar='N', type=int, help='Include powers of each confound regressor up to order N.', @@ -322,7 +322,7 @@ def _bids_filter(value, parser): ) preproc.add_argument( '--confoundderiv', - dest='confound_deriv', + dest='confoundderiv', action='store_false', help='Toggle whether derivatives will be used in confound regression.', default=True, @@ -433,7 +433,7 @@ def _bids_filter(value, parser): fixdelay = corr_fit.add_mutually_exclusive_group() fixdelay.add_argument( '--fixdelay', - dest='fixeddelayvalue', + dest='fixdelay', action='store', type=float, metavar='DELAYTIME', @@ -442,7 +442,7 @@ def _bids_filter(value, parser): ) fixdelay.add_argument( '--searchrange', - dest='lag_extrema', + dest='searchrange', action=IndicateSpecifiedAction, nargs=2, type=float, @@ -452,7 +452,7 @@ def _bids_filter(value, parser): ) corr_fit.add_argument( '--sigmalimit', - dest='widthmax', + dest='sigmalimit', action='store', type=float, metavar='SIGMALIMIT', @@ -620,7 +620,7 @@ def _bids_filter(value, parser): ) experimental.add_argument( '--autorespdelete', - dest='respdelete', + dest='autorespdelete', action='store_true', help='Attempt to detect and remove respiratory signal that strays into the LFO band.', default=False, diff --git a/src/fmripost_rapidtide/config.py b/src/fmripost_rapidtide/config.py index f6742eb..e28d949 100644 --- a/src/fmripost_rapidtide/config.py +++ b/src/fmripost_rapidtide/config.py @@ -544,21 +544,21 @@ class workflow(_Config): crosscorrelation.""" filterband = None """Filter data and regressors to specific band. Use "None" to disable filtering.""" - passvec = None + filterfreqs = None """Filter data and regressors to retain LOWERPASS to UPPERPASS.""" - stopvec = None + filterstopfreqs = None """Filter data and regressors to with stop frequencies LOWERSTOP and UPPERSTOP.""" - numestreps = None + numnull = None """Estimate significance threshold by running NREPS null correlations.""" detrendorder = None """Set order of trend removal (0 to disable).""" - gausssigma = None + spatialfilt = None """Spatially filter fMRI data prior to analysis using GAUSSSIGMA in mm.""" - confoundfilespec = None + confoundfile = None """Read additional (non-motion) confound regressors out of CONFFILE file.""" - confound_power = None + confoundpowers = None """Include powers of each confound regressor up to order N.""" - confound_deriv = None + confoundderiv = None """Toggle whether derivatives will be used in confound regression.""" globalsignalmethod = None """The method for constructing the initial global signal regressor - straight summation.""" @@ -573,11 +573,11 @@ class workflow(_Config): """Method to use for cross-correlation weighting.""" simcalcrange = None """Limit correlation calculation to data between timepoints START and END in the fmri file.""" - fixeddelayvalue = None + fixdelay = None """Don't fit the delay time - set it to DELAYTIME seconds for all voxels.""" - lag_extrema = None + searchrange = None """Limit fit to a range of lags from LAGMIN to LAGMAX.""" - widthmax = None + sigmalimit = None """Reject lag fits with linewidth wider than SIGMALIMIT Hz.""" bipolar = None """Bipolar mode - match peak correlation ignoring sign.""" @@ -605,7 +605,7 @@ class workflow(_Config): territorymap = None """This specifies a territory map. Each territory is a set of voxels with the same integral value.""" - respdelete = None + autorespdelete = None """Attempt to detect and remove respiratory signal that strays into the LFO band.""" err_on_warn = False """Cast Rapidtide warnings to errors.""" diff --git a/src/fmripost_rapidtide/data/io_spec.json b/src/fmripost_rapidtide/data/io_spec.json index a22eca0..46dfd29 100644 --- a/src/fmripost_rapidtide/data/io_spec.json +++ b/src/fmripost_rapidtide/data/io_spec.json @@ -63,17 +63,6 @@ ".nii" ] }, - "dseg": { - "datatype": "anat", - "res": "2", - "space": "MNI152NLin6Asym", - "desc": null, - "suffix": "dseg", - "extension": [ - ".nii.gz", - ".nii" - ] - }, "confounds": { "datatype": "func", "echo": null, @@ -89,6 +78,33 @@ "extension": [ ".tsv" ] + }, + "dseg_mni152nlin6asym": { + "datatype": "anat", + "task": null, + "run": null, + "res": "2", + "space": "MNI152NLin6Asym", + "desc": null, + "suffix": "dseg", + "extension": [ + ".nii.gz", + ".nii" + ] + }, + "anat_dseg": { + "datatype": "anat", + "task": null, + "run": null, + "space": null, + "res": null, + "den": null, + "desc": null, + "suffix": "dseg", + "extension": [ + ".nii.gz", + ".nii" + ] } }, "transforms": { diff --git a/src/fmripost_rapidtide/interfaces/nilearn.py b/src/fmripost_rapidtide/interfaces/nilearn.py index d58e728..1bbd928 100644 --- a/src/fmripost_rapidtide/interfaces/nilearn.py +++ b/src/fmripost_rapidtide/interfaces/nilearn.py @@ -19,7 +19,7 @@ class _MeanImageInputSpec(BaseInterfaceInputSpec): ) mask_file = File( exists=True, - mandatory=True, + mandatory=False, desc='A binary brain mask.', ) out_file = File( @@ -44,11 +44,19 @@ class MeanImage(NilearnBaseInterface, SimpleInterface): output_spec = _MeanImageOutputSpec def _run_interface(self, runtime): + import nibabel as nb from nilearn.masking import apply_mask, unmask + from nipype.interfaces.base import isdefined + + if isdefined(self.inputs.mask_file): + data = apply_mask(self.inputs.bold_file, self.inputs.mask_file) + mean_data = data.mean(axis=0) + mean_img = unmask(mean_data, self.inputs.mask_file) + else: + in_img = nb.load(self.inputs.bold_file) + mean_data = in_img.get_fdata().mean(axis=3) + mean_img = nb.Nifti1Image(mean_data, in_img.affine, in_img.header) - data = apply_mask(self.inputs.bold_file, self.inputs.mask_file) - mean_data = data.mean(axis=0) - mean_img = unmask(mean_data, self.inputs.mask_file) self._results['out_file'] = os.path.join(runtime.cwd, self.inputs.out_file) mean_img.to_filename(self._results['out_file']) diff --git a/src/fmripost_rapidtide/interfaces/rapidtide.py b/src/fmripost_rapidtide/interfaces/rapidtide.py index 213a2a7..6a3c93a 100644 --- a/src/fmripost_rapidtide/interfaces/rapidtide.py +++ b/src/fmripost_rapidtide/interfaces/rapidtide.py @@ -15,10 +15,16 @@ class _RapidtideInputSpec(CommandLineInputSpec): in_file = File( exists=True, argstr='%s', - position=-1, + position=-2, mandatory=True, desc='File to denoise', ) + outputname = traits.Str( + argstr='%s', + position=-1, + mandatory=True, + desc='Output name', + ) # Set by the workflow denoising = traits.Bool( argstr='--denoising', @@ -151,7 +157,7 @@ class _RapidtideInputSpec(CommandLineInputSpec): mandatory=False, ) timerange = traits.List( - traits.Float, + traits.Int, argstr='--timerange %s', mandatory=False, minlen=2, @@ -166,7 +172,7 @@ class _RapidtideInputSpec(CommandLineInputSpec): mandatory=False, ) simcalcrange = traits.List( - traits.Float, + traits.Int, argstr='--simcalcrange %s', mandatory=False, minlen=2, @@ -177,7 +183,7 @@ class _RapidtideInputSpec(CommandLineInputSpec): mandatory=False, ) searchrange = traits.List( - traits.Float, + traits.Int, argstr='--searchrange %s', mandatory=False, minlen=2, @@ -264,10 +270,14 @@ class Rapidtide(CommandLine): def _list_outputs(self): outputs = self._outputs().get() out_dir = os.getcwd() - outputs['delay_map'] = os.path.join(out_dir, 'desc-maxtime_map.nii.gz') + outputname = self.inputs.outputname + outputs['delay_map'] = os.path.join(out_dir, f'{outputname}_desc-maxtime_map.nii.gz') outputs['regressor_file'] = os.path.join( out_dir, - 'desc-refinedmovingregressor_timeseries.tsv.gz', + f'{outputname}_desc-refinedmovingregressor_timeseries.tsv.gz', + ) + outputs['denoised'] = os.path.join( + out_dir, + f'{outputname}_desc-lfofilterCleaned_bold.nii.gz' ) - outputs['denoised'] = os.path.join(out_dir, 'desc-lfofilterCleaned_bold.nii.gz') return outputs diff --git a/src/fmripost_rapidtide/workflows/base.py b/src/fmripost_rapidtide/workflows/base.py index 627af5e..b519b0d 100644 --- a/src/fmripost_rapidtide/workflows/base.py +++ b/src/fmripost_rapidtide/workflows/base.py @@ -383,6 +383,7 @@ def init_single_run_wf(bold_file): # Run rapidtide rapidtide_wf = init_rapidtide_wf(bold_file=bold_file, metadata=bold_metadata, mem_gb=mem_gb) rapidtide_wf.inputs.inputnode.confounds = functional_cache['confounds'] + rapidtide_wf.inputs.inputnode.dseg_std = functional_cache['dseg_mni152nlin6asym'] rapidtide_wf.inputs.inputnode.skip_vols = skip_vols mni6_buffer = pe.Node(niu.IdentityInterface(fields=['bold', 'bold_mask']), name='mni6_buffer') @@ -395,6 +396,8 @@ def init_single_run_wf(bold_file): from niworkflows.interfaces.header import ValidateImage from templateflow.api import get as get_template + raise Exception() + workflow.__desc__ += """\ Raw BOLD series were resampled to MNI152NLin6Asym:res-2, for rapidtide denoising. """ diff --git a/src/fmripost_rapidtide/workflows/rapidtide.py b/src/fmripost_rapidtide/workflows/rapidtide.py index afdacc7..2baa4e1 100644 --- a/src/fmripost_rapidtide/workflows/rapidtide.py +++ b/src/fmripost_rapidtide/workflows/rapidtide.py @@ -123,17 +123,6 @@ def init_rapidtide_wf( name='outputnode', ) - rm_non_steady_state = pe.Node( - niu.Function(function=_remove_volumes, output_names=['bold_cut']), - name='rm_nonsteady', - ) - workflow.connect([ - (inputnode, rm_non_steady_state, [ - ('skip_vols', 'skip_vols'), - ('bold_std', 'bold_file'), - ]), - ]) # fmt:skip - # Split tissue-type segmentation to get GM and WM masks split_tissues = pe.Node( SplitDseg(), @@ -142,29 +131,31 @@ def init_rapidtide_wf( workflow.connect([(inputnode, split_tissues, [('dseg_std', 'dseg')])]) # Run the Rapidtide classifier + # XXX: simcalcrange is converted to list of strings rapidtide = pe.Node( Rapidtide( + outputname='rapidtide', denoising=True, datatstep=metadata['RepetitionTime'], autosync=config.workflow.autosync, filterband=config.workflow.filterband, - filterfreqs=config.workflow.passvec or Undefined, - filterstopfreqs=config.workflow.stopvec or Undefined, - numnull=config.workflow.numestreps, + filterfreqs=config.workflow.filterfreqs or Undefined, + filterstopfreqs=config.workflow.filterstopfreqs or Undefined, + numnull=config.workflow.numnull, detrendorder=config.workflow.detrendorder, - spatialfilt=config.workflow.gausssigma, - confoundfile=config.workflow.confoundfilespec or Undefined, - confoundpowers=config.workflow.confound_power, - confoundderiv=config.workflow.confound_deriv, + spatialfilt=config.workflow.spatialfilt, + confoundfile=config.workflow.confoundfile or Undefined, + confoundpowers=config.workflow.confoundpowers, + confoundderiv=config.workflow.confoundderiv, globalsignalmethod=config.workflow.globalsignalmethod, globalpcacomponents=config.workflow.globalpcacomponents, numtozero=config.workflow.numtozero, - timerange=config.workflow.timerange, + timerange=[int(float(i)) for i in config.workflow.timerange], corrweighting=config.workflow.corrweighting, - simcalcrange=config.workflow.simcalcrange, - fixdelay=config.workflow.fixeddelayvalue or Undefined, - searchrange=config.workflow.lag_extrema, - sigmalimit=config.workflow.widthmax, + simcalcrange=[int(float(i)) for i in config.workflow.simcalcrange], + fixdelay=config.workflow.fixdelay or Undefined, + searchrange=[int(float(i)) for i in config.workflow.searchrange], + sigmalimit=config.workflow.sigmalimit, bipolar=config.workflow.bipolar, lagminthresh=config.workflow.lagminthresh, lagmaxthresh=config.workflow.lagmaxthresh, @@ -177,7 +168,7 @@ def init_rapidtide_wf( glmderivs=config.workflow.glmderivs, outputlevel=config.workflow.outputlevel, territorymap=config.workflow.territorymap or Undefined, - autorespdelete=config.workflow.respdelete, + autorespdelete=config.workflow.autorespdelete, ), name='rapidtide', mem_gb=mem_gb['resampled'], @@ -185,7 +176,7 @@ def init_rapidtide_wf( workflow.connect([ (inputnode, rapidtide, [ ('bold_std', 'in_file'), - ('brain_mask_std', 'brainmask'), + ('bold_mask_std', 'brainmask'), ('confounds', 'motionfile'), ('skip_vols', 'numskip'), ]), @@ -206,42 +197,3 @@ def init_rapidtide_wf( # Generate figures for report return workflow - - -def _remove_volumes(bold_file, skip_vols): - """Remove skip_vols from bold_file.""" - import os - - import nibabel as nb - from nipype.utils.filemanip import fname_presuffix - - if skip_vols == 0: - return bold_file - - out = fname_presuffix(bold_file, prefix='cut_', newpath=os.getcwd()) - bold_img = nb.load(bold_file) - bold_img.__class__( - bold_img.dataobj[..., skip_vols:], bold_img.affine, bold_img.header - ).to_filename(out) - return out - - -def _add_volumes(bold_file, bold_cut_file, skip_vols): - """Prepend skip_vols from bold_file onto bold_cut_file.""" - import os - - import nibabel as nb - import numpy as np - from nipype.utils.filemanip import fname_presuffix - - if skip_vols == 0: - return bold_cut_file - - bold_img = nb.load(bold_file) - bold_cut_img = nb.load(bold_cut_file) - - bold_data = np.concatenate((bold_img.dataobj[..., :skip_vols], bold_cut_img.dataobj), axis=3) - - out = fname_presuffix(bold_cut_file, prefix='addnonsteady_', newpath=os.getcwd()) - bold_img.__class__(bold_data, bold_img.affine, bold_img.header).to_filename(out) - return out