Merge pull request #1595 from FCP-INDI/feature/correlation-matrix

✨ FEAT: Correlation matrices

Author: sgiavasis

CHANGELOG.md

@@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added changelog
 - Added CHD8 mouse template (`/cpac_templates/chd8_functional_template_noise_mask_ag.nii.gz`) 
 - Added commandline flags `--T1w_label` and `--bold_label`
+- Added AFNI and Nilearn implementations of Pearson and partial correlation matrices
 
 ### Changed
 

CPAC/aroma/aroma_test.py

@@ -7,7 +7,7 @@
 file, an output_dir which can be either be mentioned or if it is set to
 None will write it in the current working directory.
 
-The argument run can either be tset to true(default) or to false. If
+The argument run can either be set to true(default) or to false. If
 set to False, it should connect to the nipype workflow and return the
 workflow object instead.
 

CPAC/connectome/connectivity_matrix.py

@@ -0,0 +1,212 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Functions for creating connectome connectivity matrices."""
+from warnings import warn
+import numpy as np
+from nilearn.connectome import ConnectivityMeasure
+from nilearn.input_data import NiftiLabelsMasker
+from nipype import logging
+from nipype.interfaces import utility as util
+from CPAC.pipeline import nipype_pipeline_engine as pe
+from CPAC.utils.interfaces.function import Function
+from CPAC.utils.interfaces.netcorr import NetCorr, strip_afni_output_header
+
+logger = logging.getLogger('nipype.workflow')
+connectome_methods = {
+    'afni': {'Pearson': '',
+             'Partial': '-part_corr'},
+    'nilearn': {'Pearson': 'correlation',
+                'Partial': 'partial correlation'}
+}
+
+
+def connectome_name(timeseries, atlas_name, tool, method):
+    """Helper function to create connectome file filename
+
+    Parameters
+    ----------
+    timeseries : str
+        path to input timeseries
+
+    atlas_name : str
+        atlas name
+
+    tool : str
+        connectome tool
+
+    method : str
+        BIDS entity value for `desc-` key
+
+    Returns
+    -------
+    str
+    """
+    method = ''.join(word.capitalize() for word in [tool, method])
+    new_filename_parts = [part for part in timeseries.split('_')[:-1][::-1] if
+                          not part.startswith('space-')]
+    atlas_index = len(new_filename_parts) - 1
+    if any(filename_part.startswith('desc-') for filename_part in
+            new_filename_parts):
+        for i, filename_part in enumerate(new_filename_parts):
+            if filename_part.startswith('desc-'):
+                new_filename_parts[-i] = f'desc-{method}'
+                atlas_index = -(i - 1)
+                break
+    new_filename_parts.insert(atlas_index, f'atlas-{atlas_name}')
+    return '_'.join([*new_filename_parts[::-1], 'connectome.tsv'])
+
+
+def get_connectome_method(method, tool):
+    """Helper function to get tool's method string
+
+    Parameters
+    ----------
+    method : str
+
+    tool : str
+
+    Returns
+    -------
+    str or NotImplemented
+
+    Examples
+    --------
+    >>> get_connectome_method('Pearson', 'AFNI')
+    ''
+    >>> get_connectome_method('Pearson', 'Nilearn')
+    'correlation'
+    >>> get_connectome_method('Spearman', 'AFNI')
+    NotImplemented
+    """
+    cm_method = connectome_methods[tool.lower()].get(method, NotImplemented)
+    if cm_method is NotImplemented:
+        warning_message = (
+            f'{method} has not yet been implemented for {tool} in C-PAC.')
+        if logger:
+            logger.warning(NotImplementedError(warning_message))
+        else:
+            warn(warning_message, category=Warning)
+    return cm_method
+
+
+def compute_connectome_nilearn(in_rois, in_file, method, atlas_name):
+    """Function to compute a connectome matrix using Nilearn
+
+    Parameters
+    ----------
+    in_rois : Niimg-like object
+        http://nilearn.github.io/manipulating_images/input_output.html#niimg-like-objects
+        Region definitions, as one image of labels.
+
+    in_file : str
+        path to timeseries image
+
+    method: str
+        'Pearson' or 'Partial'
+
+    atlas_name: str
+
+    Returns
+    -------
+    numpy.ndarray or NotImplemented
+    """
+    tool = 'Nilearn'
+    output = connectome_name(in_file, atlas_name, tool, method)
+    method = get_connectome_method(method, tool)
+    if method is NotImplemented:
+        return NotImplemented
+    masker = NiftiLabelsMasker(labels_img=in_rois,
+                               standardize=True,
+                               verbose=True)
+    timeser = masker.fit_transform(in_file)
+    correlation_measure = ConnectivityMeasure(kind=method)
+    corr_matrix = correlation_measure.fit_transform([timeser])[0]
+    np.fill_diagonal(corr_matrix, 0)
+    np.savetxt(output, corr_matrix, delimiter='\t')
+    return output
+
+
+def create_connectome_afni(name, method, pipe_num):
+    wf = pe.Workflow(name=name)
+    inputspec = pe.Node(
+        util.IdentityInterface(fields=[
+            'in_rois',  # parcellation
+            'in_file',  # timeseries,
+            'mask',
+            'method',
+            'atlas_name'
+        ]),
+        name='inputspec'
+    )
+    outputspec = pe.Node(
+        util.IdentityInterface(fields=[
+            'out_file',
+        ]),
+        name='outputspec'
+    )
+
+    timeseries_correlation = pe.Node(NetCorr(), name=name)
+    if method:
+        timeseries_correlation.inputs.part_corr = (method == 'Partial')
+
+    strip_header_node = pe.Node(Function(
+        input_names=['in_file', 'out_file'], output_names=['out_file'],
+        imports=['import subprocess'],
+        function=strip_afni_output_header),
+                                name='netcorrStripHeader'
+                                     f'{method}_{pipe_num}')
+
+    name_output_node = pe.Node(Function(input_names=['timeseries',
+                                                     'atlas_name',
+                                                     'tool',
+                                                     'method'],
+                                        output_names=['filename'],
+                                        function=connectome_name),
+                               name=f'connectomeName{method}_{pipe_num}')
+    name_output_node.inputs.tool = 'Afni'
+
+    wf.connect([
+        (inputspec, timeseries_correlation, [('in_rois', 'in_rois'),
+                                             ('in_file', 'in_file'),
+                                             ('mask', 'mask')]),
+        (inputspec, name_output_node, [('in_file', 'timeseries'),
+                                       ('atlas_name', 'atlas_name'),
+                                       ('method', 'method')]),
+        (timeseries_correlation, strip_header_node, [
+            ('out_corr_matrix', 'in_file')]),
+        (name_output_node, strip_header_node, [('filename', 'out_file')]),
+        (strip_header_node, outputspec, [('out_file', 'out_file')])])
+    return wf
+
+
+def create_connectome_nilearn(name='connectomeNilearn'):
+    wf = pe.Workflow(name=name)
+    inputspec = pe.Node(
+        util.IdentityInterface(fields=[
+            'in_rois',  # parcellation
+            'in_file',  # timeseries
+            'method',
+            'atlas_name'
+        ]),
+        name='inputspec'
+    )
+    outputspec = pe.Node(
+        util.IdentityInterface(fields=[
+            'out_file',
+        ]),
+        name='outputspec'
+    )
+    node = pe.Node(Function(input_names=['in_rois', 'in_file', 'method',
+                                         'atlas_name'],
+                            output_names=['out_file'],
+                            function=compute_connectome_nilearn,
+                            as_module=True),
+                   name='connectome')
+    wf.connect([
+        (inputspec, node, [('in_rois', 'in_rois'),
+                           ('in_file', 'in_file'),
+                           ('method', 'method'),
+                           ('atlas_name', 'atlas_name')]),
+        (node, outputspec, [('out_file', 'out_file')]),
+    ])
+    return wf

CPAC/connectome/pipeline.py

@@ -0,0 +1,8 @@
+"""Tests for CPAC.connectome"""
+import os
+import random
+import re
+import numpy as np
+import pytest
+from CPAC.pipeline.schema import valid_options
+

CPAC/connectome/test/test_connectome.py

@@ -25,7 +25,7 @@ def run_warp_nipype(inputs, output_dir=None, run=True):
     #  warp_nipe file, an output_dir which can be either be mentioned
     #  or if it is set to none will write it in the current working
     #  directory.
-    #  The argument run can either be tset to true(default) or to
+    #  The argument run can either be set to true(default) or to
     #  false. If set to false, it should connect to the nipype workflow
     #  and return the workflow object instead
     #  What all should it return?: if the run had been set to true, it

CPAC/distortion_correction/tests/test_distortion_correction.py

@@ -984,7 +984,36 @@ def build_segmentation_stack(rpool, cfg, pipeline_blocks=None):
     return pipeline_blocks
 
 
+def list_blocks(pipeline_blocks, indent=None):
+    """Function to list node blocks line by line
+
+    Parameters
+    ----------
+    pipeline_blocks : list or tuple
+
+    indent : int or None
+       number of spaces after a tab indent
+
+    Returns
+    -------
+    str
+    """
+    blockstring = yaml.dump([
+        getattr(block, '__name__', getattr(block, 'name', yaml.safe_load(
+            list_blocks(list(block))) if
+            isinstance(block, (tuple, list, set)) else str(block))
+        ) for block in pipeline_blocks
+    ])
+    if isinstance(indent, int):
+        blockstring = '\n'.join([
+            '\t' + ' ' * indent + line for line in blockstring.split('\n')])
+    return blockstring
+
+
 def connect_pipeline(wf, cfg, rpool, pipeline_blocks):
+    logger.info('\n'.join([
+        'Connecting pipeline blocks:',
+        list_blocks(pipeline_blocks, indent=1)]))
 
     for block in pipeline_blocks:
         try:

CPAC/pipeline/cpac_pipeline.py

@@ -42,8 +42,13 @@
         'template': ['EPI_Template', 'T1_Template'],
     },
     'timeseries': {
-        'roi_paths': {'Avg', 'Voxel', 'SpatialReg', 'PearsonCorr',
-                      'PartialCorr'},
+        'roi_paths': {'Avg', 'Voxel', 'SpatialReg'},
+    },
+    'connectivity_matrix': {
+        'using': {'AFNI', 'Nilearn'},
+        'measure': {'Pearson', 'Partial', 'Spearman', 'MGC',
+                    # 'TangentEmbed'  # "Skip tangent embedding for now"
+                    },
     },
     'Regressors': {
         'CompCor': {
@@ -965,8 +970,11 @@ def _changes_1_8_0_to_1_8_1(config_dict):
                 'tse_roi_paths', valid_options['timeseries']['roi_paths'])
         ),
         'realignment': In({'ROI_to_func', 'func_to_ROI'}),
+        'connectivity_matrix': {
+            option: Maybe([In(valid_options['connectivity_matrix'][option])])
+            for option in ['using', 'measure']
+        },
     },
-
     'seed_based_correlation_analysis': {
         'run': bool,
         Optional('roi_paths_fully_specified'): bool,

CPAC/pipeline/schema.py

@@ -206,7 +206,7 @@ timeseries_extraction:
   # Enter paths to region-of-interest (ROI) NIFTI files (.nii or .nii.gz) to be used for time-series extraction, and then select which types of analyses to run.
   # Denote which analyses to run for each ROI path by listing the names below. For example, if you wish to run Avg and SpatialReg, you would enter: '/path/to/ROI.nii.gz': Avg, SpatialReg
   # available analyses:
-  #   /path/to/atlas.nii.gz: Avg, Voxel, SpatialReg, PearsonCorr, PartialCorr
+  #   /path/to/atlas.nii.gz: Avg, Voxel, SpatialReg
   tse_roi_paths: 
     s3://fcp-indi/resources/cpac/resources/rois_2mm.nii.gz: Avg, Voxel, SpatialReg
 

CPAC/resources/configs/pipeline_config_benchmark-ANTS.yml

@@ -225,7 +225,7 @@ timeseries_extraction:
   # Enter paths to region-of-interest (ROI) NIFTI files (.nii or .nii.gz) to be used for time-series extraction, and then select which types of analyses to run.
   # Denote which analyses to run for each ROI path by listing the names below. For example, if you wish to run Avg and SpatialReg, you would enter: '/path/to/ROI.nii.gz': Avg, SpatialReg
   # available analyses:
-  #   /path/to/atlas.nii.gz: Avg, Voxel, SpatialReg, PearsonCorr, PartialCorr
+  #   /path/to/atlas.nii.gz: Avg, Voxel, SpatialReg
   tse_roi_paths: 
     s3://fcp-indi/resources/cpac/resources/rois_2mm.nii.gz: Avg, Voxel, SpatialReg