Add orthogonalized motion components to confounds file (#56)

* Add orthogonalized components to confounds.

* Replace unstable with main.

.circleci/config.yml

@@ -274,8 +274,8 @@ jobs:
           command: |
             if [[ -n "$DOCKER_PASS" ]]; then
               docker login -u $DOCKER_USER -p $DOCKER_PASS
-              docker tag nipreps/fmripost_aroma nipreps/fmripost_aroma:unstable
-              docker push nipreps/fmripost_aroma:unstable
+              docker tag nipreps/fmripost_aroma nipreps/fmripost_aroma:main
+              docker push nipreps/fmripost_aroma:main
               if [[ -n "$CIRCLE_TAG" ]]; then
                 docker push nipreps/fmripost_aroma:latest
                 docker tag nipreps/fmripost_aroma nipreps/fmripost_aroma:$CIRCLE_TAG

README.md

@@ -20,7 +20,7 @@ to MNI152NLin6Asym at 2 mm3 resolution.
 ## Installation
 
 ```console
-docker pull nipreps/fmripost-aroma:unstable
+docker pull nipreps/fmripost-aroma:main
 ```
 
 ## License

docs/installation.rst

@@ -6,7 +6,7 @@ Installation
 *fMRIPost-AROMA* should be installed using container technologies.
 
 .. code-block:: bash
-  docker pull nipreps/fmripost-aroma:unstable
+  docker pull nipreps/fmripost-aroma:main
 
 
 Containerized execution (Docker and Singularity)

docs/outputs.rst

@@ -24,8 +24,7 @@ upcoming `BEP 011`_ and `BEP 012`_).
     Denoised fMRI data in the requested output spaces and resolutions.
 
 4.  **Confounds**:
-    Time series of ICA components,
-    as well as labels indicating if the components are accepted or rejected.
+    Time series of ICA components classified as noise.
 
 
 Layout
@@ -107,19 +106,23 @@ Confounds_ are saved as a :abbr:`TSV (tab-separated value)` file::
     func/
       sub-<label>_[specifiers]_desc-aroma_metrics.tsv
       sub-<label>_[specifiers]_desc-aroma_metrics.json
-      sub-<label>_[specifiers]_desc-melodic_timeseries.tsv
-      sub-<label>_[specifiers]_desc-melodic_timeseries.json
+      sub-<label>_[specifiers]_desc-aroma_timeseries.tsv
+      sub-<label>_[specifiers]_desc-aroma_timeseries.json
 
 Confounds
 ---------
 
 *fMRIPost-AROMA* outputs a set of confounds that can be used to denoise the data.
-These are stored in a TSV file (``desc-melodic_timeseries.tsv``) and a JSON file
-(``desc-melodic_timeseries.json``) that contains metadata about the confounds.
+These are stored in a TSV file (``desc-aroma_timeseries.tsv``) and a JSON file
+(``desc-aroma_timeseries.json``) that contains metadata about the confounds.
 
 The confounds generated by *fMRIPost-AROMA* are ICA component time series
 classified as "rejected" by ICA-AROMA.
 
+Columns starting with ``aroma_motion_`` are the raw noise ICA component time series.
+Columns starting with ``aroma_orth_motion_`` are the noise ICA component time series,
+after z-scoring and orthogonalization with respect to the signal ICA component time series.
+
 Confounds and "carpet"-plot on the visual reports
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

src/fmripost_aroma/interfaces/confounds.py

@@ -68,6 +68,9 @@ def _get_ica_confounds(mixing, aroma_features, skip_vols, newpath=None):
     motion_ics = aroma_features_df.loc[
         aroma_features_df['classification'] == 'rejected'
     ].index.values
+    signal_ics = aroma_features_df.loc[
+        aroma_features_df['classification'] != 'rejected'
+    ].index.values
     mixing_arr = np.loadtxt(mixing, ndmin=2)
 
     # Prepare output paths
@@ -96,11 +99,29 @@ def _get_ica_confounds(mixing, aroma_features, skip_vols, newpath=None):
     # Select the mixing matrix rows corresponding to the motion ICs
     aggr_mixing_arr = mixing_arr[motion_ics, :].T
 
+    signal_mixing_arr = mixing_arr[signal_ics, :].T
+    orthaggr_mixing_arr = aggr_mixing_arr.copy()
+    orthaggr_mixing_arr = aggr_mixing_arr - np.dot(
+        np.dot(np.linalg.pinv(good_ic_arr), good_ic_arr), aggr_mixing_arr
+    )
+    # Regress the good components out of the bad time series to get "pure evil" regressors
+    aggr_mixing_arr_z = stats.zscore(aggr_mixing_arr, axis=0)
+    signal_mixing_arr_z = stats.zscore(signal_mixing_arr, axis=0)
+    betas = np.linalg.lstsq(signal_mixing_arr_z, aggr_mixing_arr_z, rcond=None)[0]
+    pred_bad_timeseries = np.dot(signal_mixing_arr_z, betas)
+    orthaggr_mixing_arr = aggr_mixing_arr_z - pred_bad_timeseries
+
     # add one to motion_ic_indices to match melodic report.
-    pd.DataFrame(
+    aggr_confounds_df = pd.DataFrame(
         aggr_mixing_arr,
         columns=[f'aroma_motion_{x + 1:02d}' for x in motion_ics],
-    ).to_csv(aroma_confounds, sep='\t', index=None)
+    )
+    orthaggr_confounds_df = pd.DataFrame(
+        orthaggr_mixing_arr,
+        columns=[f'aroma_orth_motion_{x + 1:02d}' for x in motion_ics],
+    )
+    confounds_df = pd.concat([aggr_confounds_df, orthaggr_confounds_df], axis=1)
+    confounds_df.to_csv(aroma_confounds, sep='\t', index=None)
 
     return aroma_confounds, mixing_out
 

src/fmripost_aroma/tests/data/test_ds005115_deriv_and_raw_outputs.txt

@@ -4,7 +4,7 @@ sub-01/ses-01
 sub-01/ses-01/func
 sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-aroma_metrics.json
 sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-aroma_metrics.tsv
-sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-melodic_timeseries.tsv
+sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-aroma_timeseries.tsv
 sub-01/ses-01/func/sub-01_ses-01_task-rest_space-MNI152NLin6Asym_res-2_desc-aggrDenoised_bold.nii.gz
 sub-01/ses-01/func/sub-01_ses-01_task-rest_space-MNI152NLin6Asym_res-2_desc-melodic_components.nii.gz
 sub-01/ses-01/func/sub-01_ses-01_task-rest_space-MNI152NLin6Asym_res-2_desc-melodic_mixing.tsv

src/fmripost_aroma/tests/data/test_ds005115_deriv_only_outputs.txt

@@ -4,7 +4,7 @@ sub-01/ses-01
 sub-01/ses-01/func
 sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-aroma_metrics.json
 sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-aroma_metrics.tsv
-sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-melodic_timeseries.tsv
+sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-aroma_timeseries.tsv
 sub-01/ses-01/func/sub-01_ses-01_task-rest_space-MNI152NLin6Asym_res-2_desc-aggrDenoised_bold.nii.gz
 sub-01/ses-01/func/sub-01_ses-01_task-rest_space-MNI152NLin6Asym_res-2_desc-melodic_components.nii.gz
 sub-01/ses-01/func/sub-01_ses-01_task-rest_space-MNI152NLin6Asym_res-2_desc-melodic_mixing.tsv

src/fmripost_aroma/tests/data/test_ds005115_resampling_and_raw_outputs.txt

@@ -4,7 +4,7 @@ sub-01/ses-01
 sub-01/ses-01/func
 sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-aroma_metrics.json
 sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-aroma_metrics.tsv
-sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-melodic_timeseries.tsv
+sub-01/ses-01/func/sub-01_ses-01_task-rest_desc-aroma_timeseries.tsv
 sub-01/ses-01/func/sub-01_ses-01_task-rest_space-MNI152NLin6Asym_res-2_desc-aggrDenoised_bold.nii.gz
 sub-01/ses-01/func/sub-01_ses-01_task-rest_space-MNI152NLin6Asym_res-2_desc-melodic_components.nii.gz
 sub-01/ses-01/func/sub-01_ses-01_task-rest_space-MNI152NLin6Asym_res-2_desc-melodic_mixing.tsv

src/fmripost_aroma/tests/run_local_tests.py

@@ -73,7 +73,7 @@ def run_tests(test_regex, test_mark):
     run_str = 'docker run --rm -ti '
     run_str += f'-v {local_patch}:{mounted_code} '
     run_str += '--entrypoint pytest '
-    run_str += 'nipreps/fmripost_aroma:unstable '
+    run_str += 'nipreps/fmripost_aroma:main '
     run_str += (
         f'{mounted_code}/fmripost_aroma '
         f'--data_dir={mounted_code}/fmripost_aroma/tests/test_data '

src/fmripost_aroma/workflows/aroma.py

@@ -352,7 +352,7 @@ def init_ica_aroma_wf(
             base_directory=config.execution.output_dir,
             source_file=bold_file,
             datatype='func',
-            desc='melodic',
+            desc='aroma',
             suffix='timeseries',
             extension='tsv',
             dismiss_entities=('echo', 'den', 'res'),