add headers to outputs of compcor, framewise displacement + test fix

Author: Shoshana Berleant

nipype/algorithms/confounds.py

@@ -255,7 +255,7 @@ def _run_interface(self, runtime):
             'out_file': op.abspath(self.inputs.out_file),
             'fd_average': float(fd_res.mean())
         }
-        np.savetxt(self.inputs.out_file, fd_res)
+        np.savetxt(self.inputs.out_file, fd_res, header='framewise_displacement')
 
         if self.inputs.save_plot:
             tr = None
@@ -291,6 +291,8 @@ class CompCorInputSpec(BaseInterfaceInputSpec):
                                    'pre-component extraction')
     regress_poly_degree = traits.Range(low=1, default=1, usedefault=True,
                                        desc='the degree polynomial to use')
+    header = traits.Str(desc='the desired header for the output tsv file (one column).'
+                        'If undefined, will default to "CompCor"')
 
 class CompCorOutputSpec(TraitedSpec):
     components_file = File(exists=True,
@@ -359,7 +361,9 @@ def _run_interface(self, runtime):
         u, _, _ = linalg.svd(M, full_matrices=False)
         components = u[:, :self.inputs.num_components]
         components_file = os.path.join(os.getcwd(), self.inputs.components_file)
-        np.savetxt(components_file, components, fmt=b"%.10f")
+
+        self._set_header()
+        np.savetxt(components_file, components, fmt=b"%.10f", header=self._make_headers(components.shape[1]))
         return runtime
 
     def _list_outputs(self):
@@ -374,6 +378,26 @@ def _compute_tSTD(self, M, x):
         stdM[np.isnan(stdM)] = x
         return stdM
 
+    def _set_header(self, header='CompCor'):
+        self.inputs.header = self.inputs.header if isdefined(self.inputs.header) else header
+
+    def _make_headers(self, num_col):
+        headers = []
+        for i in range(num_col):
+            headers.append(self.inputs.header + str(i))
+        return '\t'.join(headers)
+
+
+class ACompCor(CompCor):
+    ''' Anatomical compcor; for input/output, see CompCor.
+    If the mask provided is an anatomical mask, CompCor == ACompCor '''
+
+    def __init__(self, *args, **kwargs):
+        ''' exactly the same as compcor except the header '''
+        super(ACompCor, self).__init__(*args, **kwargs)
+        self._set_header('aCompCor')
+
+
 class TCompCorInputSpec(CompCorInputSpec):
     # and all the fields in CompCorInputSpec
     percentile_threshold = traits.Range(low=0., high=1., value=.02,
@@ -439,12 +463,11 @@ def _run_interface(self, runtime):
         IFLOG.debug('tCompcor computed and saved mask of shape {} to mask_file {}'
                    .format(mask.shape, mask_file))
         self.inputs.mask_file = mask_file
+        self._set_header('tCompCor')
 
         super(TCompCor, self)._run_interface(runtime)
         return runtime
 
-ACompCor = CompCor
-
 class TSNRInputSpec(BaseInterfaceInputSpec):
     in_file = InputMultiPath(File(exists=True), mandatory=True,
                              desc='realigned 4D file or a list of 3D files')

nipype/algorithms/tests/test_compcor.py

@@ -38,25 +38,20 @@ def test_compcor(self):
                                ['0.4206466244', '-0.3361270124'],
                                ['-0.1246655485', '-0.1235705610']]
 
-        ccresult = self.run_cc(CompCor(realigned_file=self.realigned_file,
-                                       mask_file=self.mask_file),
+        self.run_cc(CompCor(realigned_file=self.realigned_file, mask_file=self.mask_file),
                                expected_components)
 
-        accresult = self.run_cc(ACompCor(realigned_file=self.realigned_file,
-                                         mask_file=self.mask_file,
-                                         components_file='acc_components_file'),
-                                expected_components)
-
-        assert_equal(os.path.getsize(ccresult.outputs.components_file),
-                     os.path.getsize(accresult.outputs.components_file))
+        self.run_cc(ACompCor(realigned_file=self.realigned_file, mask_file=self.mask_file,
+                             components_file='acc_components_file'),
+                    expected_components, 'aCompCor')
 
     def test_tcompcor(self):
         ccinterface = TCompCor(realigned_file=self.realigned_file, percentile_threshold=0.75)
         self.run_cc(ccinterface, [['-0.1114536190', '-0.4632908609'],
                                   ['0.4566907310', '0.6983205193'],
                                   ['-0.7132557407', '0.1340170559'],
                                   ['0.5022537643', '-0.5098322262'],
-                                  ['-0.1342351356', '0.1407855119']])
+                                  ['-0.1342351356', '0.1407855119']], 'tCompCor')
 
     def test_tcompcor_no_percentile(self):
         ccinterface = TCompCor(realigned_file=self.realigned_file)
@@ -96,7 +91,7 @@ def test_tcompcor_bad_input_dim(self):
         interface = TCompCor(realigned_file=data_file)
         self.assertRaisesRegexp(ValueError, '4-D', interface.run)
 
-    def run_cc(self, ccinterface, expected_components):
+    def run_cc(self, ccinterface, expected_components, expected_header='CompCor'):
         # run
         ccresult = ccinterface.run()
 
@@ -108,12 +103,19 @@ def run_cc(self, ccinterface, expected_components):
         assert_equal(ccinterface.inputs.num_components, 6)
 
         with open(ccresult.outputs.components_file, 'r') as components_file:
+            expected_n_components = min(ccinterface.inputs.num_components, self.fake_data.shape[3])
+
             components_data = [line.split() for line in components_file]
-            num_got_components = len(components_data)
-            assert_true(num_got_components == ccinterface.inputs.num_components
-                        or num_got_components == self.fake_data.shape[3])
-            first_two = [row[:2] for row in components_data]
-            assert_equal(first_two, expected_components)
+
+            header = components_data.pop(0)[1:] # the first item will be '#', we can throw it out
+            assert_equal(header, [expected_header + str(i) for i in range(expected_n_components)])
+
+            num_got_timepoints = len(components_data)
+            assert_equal(num_got_timepoints, self.fake_data.shape[3])
+            for index, timepoint in enumerate(components_data):
+                assert_true(len(timepoint) == ccinterface.inputs.num_components
+                            or len(timepoint) == self.fake_data.shape[3])
+                assert_equal(timepoint[:2], expected_components[index])
         return ccresult
 
     def tearDown(self):

nipype/algorithms/tests/test_confounds.py

@@ -4,7 +4,7 @@
 from tempfile import mkdtemp
 from shutil import rmtree
 
-from nipype.testing import (assert_equal, example_data, skipif, assert_true)
+from nipype.testing import (assert_equal, example_data, skipif, assert_true, assert_in)
 from nipype.algorithms.confounds import FramewiseDisplacement, ComputeDVARS
 import numpy as np
 
@@ -24,8 +24,14 @@ def test_fd():
                                           out_file=tempdir + '/fd.txt')
     res = fdisplacement.run()
 
+    with open(res.outputs.out_file) as all_lines:
+        for line in all_lines:
+            yield assert_in, 'framewise_displacement', line
+            break
+
     yield assert_true, np.allclose(ground_truth, np.loadtxt(res.outputs.out_file), atol=.16)
     yield assert_true, np.abs(ground_truth.mean() - res.outputs.fd_average) < 1e-2
+
     rmtree(tempdir)
 
 @skipif(nonitime)
@@ -40,3 +46,5 @@ def test_dvars():
 
     dv1 = np.loadtxt(res.outputs.out_std)
     yield assert_equal, (np.abs(dv1 - ground_truth).sum()/ len(dv1)) < 0.05, True
+
+    rmtree(tempdir)