Keep working.

src/fmripost_aroma/config.py

@@ -396,6 +396,8 @@ class execution(_Config):
     """Debug mode(s)."""
     fmripost_aroma_dir = None
     """Root of fMRIPost-AROMA BIDS Derivatives dataset. Depends on output_layout."""
+    fs_license_file = _fs_license
+    """An existing file containing a FreeSurfer license."""
     layout = None
     """A :py:class:`~bids.layout.BIDSLayout` object, see :py:func:`init`."""
     log_dir = None
@@ -441,6 +443,7 @@ class execution(_Config):
         'derivatives',
         'bids_database_dir',
         'fmripost_aroma_dir',
+        'fs_license_file',
         'layout',
         'log_dir',
         'output_dir',

src/fmripost_aroma/data/tests/config.toml

@@ -30,8 +30,8 @@ write_graph = false
 
 [workflow]
 anat_only = false
-aroma_err_on_warn = false
-aroma_melodic_dim = -200
+err_on_warn = false
+melodic_dim = -200
 bold2anat_dof = 6
 fmap_bspline = false
 force_syn = false

src/fmripost_aroma/interfaces/aroma.py

@@ -21,12 +21,12 @@ class _AROMAClassifierInputSpec(BaseInterfaceInputSpec):
     motpars = File(exists=True, desc='motion parameters or general confounds file')
     mixing = File(exists=True, desc='mixing matrix')
     component_maps = File(exists=True, desc='thresholded z-statistic component maps')
+    component_stats = File(exists=True, desc='melodic component statistics file')
     TR = traits.Float(desc='repetition time in seconds')
 
 
 class _AROMAClassifierOutputSpec(TraitedSpec):
     aroma_features = File(exists=True, desc='output confounds file extracted from ICA-AROMA')
-    aroma_noise_ics = File(exists=True, desc='ICA-AROMA noise components')
     aroma_metadata = traits.Dict(desc='metadata for the ICA-AROMA confounds')
 
 
@@ -77,19 +77,24 @@ def _run_interface(self, runtime):
         features_df = pd.concat([features_df, clf_df], axis=1)
         metric_metadata.update(clf_metadata)
 
-        # Split out the motion components
-        motion_components = features_df.loc[
-            features_df['classification'] == 'rejected',
-            'classification',
-        ].index
-        motion_components = motion_components.values
-        motion_components_file = os.path.abspath('aroma_motion_components.txt')
-        np.savetxt(motion_components_file, motion_components, fmt='%s')
+        # Add MELODIC component statistics to the AROMA features
+        component_stats = pd.read_csv(component_stats, header=None, sep='  ')[[0, 1]] / 100
+        component_stats.columns = ['model_variance_explained', 'total_variance_explained']
+        features_df = pd.concat([features_df, component_stats], axis=1)
+        metric_metadata.update(
+            {
+                'model_variance_explained': {
+                    'Description': 'Model variance explained by the component',
+                },
+                'total_variance_explained': {
+                    'Description': 'Total variance explained by the component',
+                },
+            },
+        )
 
         features_file = os.path.abspath('aroma_features.tsv')
         features_df.to_csv(features_file, sep='\t', index=False)
 
         self._results['aroma_features'] = features_file
         self._results['aroma_metadata'] = metric_metadata
-        self._results['aroma_noise_ics'] = motion_components_file
         return runtime

src/fmripost_aroma/interfaces/confounds.py

@@ -1,14 +1,11 @@
 """Handling confounds."""
 
 import os
-import re
-import shutil
 
 import numpy as np
 import pandas as pd
 from nipype.interfaces.base import (
     BaseInterfaceInputSpec,
-    Directory,
     File,
     SimpleInterface,
     TraitedSpec,
@@ -19,10 +16,8 @@
 
 
 class _ICAConfoundsInputSpec(BaseInterfaceInputSpec):
-    in_directory = Directory(
-        mandatory=True,
-        desc='directory where ICA derivatives are found',
-    )
+    mixing = File(exists=True, desc='melodic mixing matrix')
+    aroma_features = File(exists=True, desc='output confounds file extracted from ICA-AROMA')
     skip_vols = traits.Int(desc='number of non steady state volumes identified')
     err_on_aroma_warn = traits.Bool(False, usedefault=True, desc='raise error if aroma fails')
 
@@ -32,8 +27,7 @@ class _ICAConfoundsOutputSpec(TraitedSpec):
         None, File(exists=True, desc='output confounds file extracted from ICA-AROMA')
     )
     aroma_noise_ics = File(exists=True, desc='ICA-AROMA noise components')
-    melodic_mix = File(exists=True, desc='melodic mix file')
-    aroma_metadata = File(exists=True, desc='tabulated ICA-AROMA metadata')
+    mixing = File(exists=True, desc='melodic mix file with skip_vols added back')
 
 
 class ICAConfounds(SimpleInterface):
@@ -43,89 +37,74 @@ class ICAConfounds(SimpleInterface):
     output_spec = _ICAConfoundsOutputSpec
 
     def _run_interface(self, runtime):
-        (aroma_confounds, motion_ics_out, melodic_mix_out, aroma_metadata) = _get_ica_confounds(
-            self.inputs.in_directory, self.inputs.skip_vols, newpath=runtime.cwd
+        aroma_confounds, motion_ics_out, mixing = _get_ica_confounds(
+            mixing=self.inputs.mixing,
+            aroma_features=self.inputs.aroma_features,
+            skip_vols=self.inputs.skip_vols,
+            newpath=runtime.cwd,
         )
 
         if self.inputs.err_on_aroma_warn and aroma_confounds is None:
             raise RuntimeError('ICA-AROMA failed')
 
-        aroma_confounds = self._results['aroma_confounds'] = aroma_confounds
-
+        self._results['aroma_confounds'] = aroma_confounds
         self._results['aroma_noise_ics'] = motion_ics_out
-        self._results['melodic_mix'] = melodic_mix_out
-        self._results['aroma_metadata'] = aroma_metadata
+        self._results['mixing'] = mixing
         return runtime
 
 
-def _get_ica_confounds(ica_out_dir, skip_vols, newpath=None):
-    """Extract confounds from ICA-AROMA result directory."""
+def _get_ica_confounds(mixing, aroma_features, skip_vols, newpath=None):
+    """Extract confounds from ICA-AROMA result directory.
+
+    This function does the following:
+
+    1. Add the number of non steady state volumes to the melodic mix file.
+    2. Select motion ICs from the mixing matrix for new "confounds" file.
+    """
     if newpath is None:
         newpath = os.getcwd()
 
-    # load the txt files from ICA-AROMA
-    melodic_mix = os.path.join(ica_out_dir, 'melodic.ica/melodic_mix')
-    motion_ics = os.path.join(ica_out_dir, 'classified_motion_ICs.txt')
-    aroma_metadata = os.path.join(ica_out_dir, 'classification_overview.txt')
-    aroma_icstats = os.path.join(ica_out_dir, 'melodic.ica/melodic_ICstats')
+    # Load input files
+    aroma_features_df = pd.read_table(aroma_features)
+    motion_ics = aroma_features_df.loc[
+        aroma_features_df['classification'] == 'rejected'
+    ].index.values
+    mixing_arr = np.loadtxt(mixing, ndmin=2)
 
-    # Change names of motion_ics and melodic_mix for output
-    melodic_mix_out = os.path.join(newpath, 'MELODICmix.tsv')
+    # Prepare output paths
+    mixing_out = os.path.join(newpath, 'mixing.tsv')
     motion_ics_out = os.path.join(newpath, 'AROMAnoiseICs.csv')
-    aroma_metadata_out = os.path.join(newpath, 'classification_overview.tsv')
-
-    # copy metion_ics file to derivatives name
-    shutil.copyfile(motion_ics, motion_ics_out)
-
-    # -1 since python lists start at index 0
-    motion_ic_indices = np.loadtxt(motion_ics, dtype=int, delimiter=',', ndmin=1) - 1
-    melodic_mix_arr = np.loadtxt(melodic_mix, ndmin=2)
+    aroma_confounds = os.path.join(newpath, 'AROMAAggrCompAROMAConfounds.tsv')
 
-    # pad melodic_mix_arr with rows of zeros corresponding to number non steadystate volumes
+    # pad mixing_arr with rows of zeros corresponding to number non steady-state volumes
     if skip_vols > 0:
-        zeros = np.zeros([skip_vols, melodic_mix_arr.shape[1]])
-        melodic_mix_arr = np.vstack([zeros, melodic_mix_arr])
-
-    # save melodic_mix_arr
-    np.savetxt(melodic_mix_out, melodic_mix_arr, delimiter='\t')
-
-    # process the metadata so that the IC column entries match the BIDS name of
-    # the regressor
-    aroma_metadata = pd.read_csv(aroma_metadata, sep='\t')
-    aroma_metadata['IC'] = [f'aroma_motion_{name}' for name in aroma_metadata['IC']]
-    aroma_metadata.columns = [re.sub(r'[ |\-|\/]', '_', c) for c in aroma_metadata.columns]
+        zeros = np.zeros([skip_vols, mixing_arr.shape[1]])
+        mixing_arr = np.vstack([zeros, mixing_arr])
 
-    # Add variance statistics to metadata
-    aroma_icstats = pd.read_csv(aroma_icstats, header=None, sep='  ')[[0, 1]] / 100
-    aroma_icstats.columns = ['model_variance_explained', 'total_variance_explained']
-    aroma_metadata = pd.concat([aroma_metadata, aroma_icstats], axis=1)
-
-    aroma_metadata.to_csv(aroma_metadata_out, sep='\t', index=False)
+    # save mixing_arr
+    np.savetxt(mixing_out, mixing_arr, delimiter='\t')
 
     # Return dummy list of ones if no noise components were found
-    if motion_ic_indices.size == 0:
+    if motion_ics.size == 0:
         config.loggers.interfaces.warning('No noise components were classified')
-        return None, motion_ics_out, melodic_mix_out, aroma_metadata_out
-
-    # the "good" ics, (e.g., not motion related)
-    good_ic_arr = np.delete(melodic_mix_arr, motion_ic_indices, 1).T
+        return None, motion_ics_out, mixing_out
 
     # return dummy lists of zeros if no signal components were found
+    good_ic_arr = np.delete(mixing_arr, motion_ics, 1).T
     if good_ic_arr.size == 0:
         config.loggers.interfaces.warning('No signal components were classified')
-        return None, motion_ics_out, melodic_mix_out, aroma_metadata_out
+        return None, motion_ics_out, mixing_out
 
-    # transpose melodic_mix_arr so x refers to the correct dimension
-    aggr_confounds = np.asarray([melodic_mix_arr.T[x] for x in motion_ic_indices])
+    # Select the mixing matrix rows corresponding to the motion ICs
+    aggr_mixing_arr = mixing_arr[motion_ics, :].T
 
     # add one to motion_ic_indices to match melodic report.
-    aroma_confounds = os.path.join(newpath, 'AROMAAggrCompAROMAConfounds.tsv')
     pd.DataFrame(
-        aggr_confounds.T,
-        columns=[f'aroma_motion_{x + 1:02d}' for x in motion_ic_indices],
+        aggr_mixing_arr,
+        columns=[f'aroma_motion_{x + 1:02d}' for x in motion_ics],
     ).to_csv(aroma_confounds, sep='\t', index=None)
 
-    return aroma_confounds, motion_ics_out, melodic_mix_out, aroma_metadata_out
+    return aroma_confounds, motion_ics_out, mixing_out
 
 
 class _ICADenoiseInputSpec(BaseInterfaceInputSpec):