From 7e2837be0c06e22a263ae5b6709649602f12aba6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Wed, 10 Apr 2024 12:38:30 +0200 Subject: [PATCH 01/12] Starting DC61 repro [skip ci] --- narps_open/pipelines/__init__.py | 2 +- narps_open/pipelines/team_DC61.py | 690 ++++++++++++++++++ tests/pipelines/test_team_DC61.py | 128 ++++ .../pipelines/team_DC61/confounds.tsv | 3 + 4 files changed, 822 insertions(+), 1 deletion(-) create mode 100644 narps_open/pipelines/team_DC61.py create mode 100644 tests/pipelines/test_team_DC61.py create mode 100644 tests/test_data/pipelines/team_DC61/confounds.tsv diff --git a/narps_open/pipelines/__init__.py b/narps_open/pipelines/__init__.py index b8ca93d4..dd6f5a38 100644 --- a/narps_open/pipelines/__init__.py +++ b/narps_open/pipelines/__init__.py @@ -44,7 +44,7 @@ 'B5I6': None, 'C22U': None, 'C88N': 'PipelineTeamC88N', - 'DC61': None, + 'DC61': 'PipelineTeamDC61', 'E3B6': None, 'E6R3': None, 'I07H': None, diff --git a/narps_open/pipelines/team_DC61.py b/narps_open/pipelines/team_DC61.py new file mode 100644 index 00000000..ba04a626 --- /dev/null +++ b/narps_open/pipelines/team_DC61.py @@ -0,0 +1,690 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Write the work of NARPS' team DC61 using Nipype """ + +from os.path import join +from itertools import product + +from nipype import Workflow, Node, MapNode +from nipype.interfaces.utility import IdentityInterface, Function +from nipype.interfaces.io import SelectFiles, DataSink +from nipype.interfaces.spm import ( + Smooth, + OneSampleTTestDesign, EstimateModel, EstimateContrast, + Level1Design, TwoSampleTTestDesign, Threshold + ) +from nipype.algorithms.modelgen import SpecifySPMModel +from nipype.algorithms.misc import Gunzip + +from narps_open.pipelines import Pipeline +from narps_open.data.task import TaskInformation +from narps_open.data.participants import get_group +from narps_open.core.interfaces import InterfaceFactory +from narps_open.core.common import list_intersection, elements_in_string, clean_list +from narps_open.utils.configuration import Configuration + +class PipelineTeamDC61(Pipeline): + """ A class that defines the pipeline of team DC61. """ + + def __init__(self): + super().__init__() + self.fwhm = 5.0 + self.team_id = 'DC61' + self.contrast_list = ['0001', '0002'] + self.subject_level_contrasts = [ + ('effect_of_gain', 'T', + [f'gamble_run{r}xgain_param^1' for r in range(1, len(self.run_list) + 1)], + [0.25]*len(self.run_list)), + ('effect_of_loss', 'T', + [f'gamble_run{r}xloss_param^1' for r in range(1, len(self.run_list) + 1)], + [0.25]*len(self.run_list)) + ] + + def get_preprocessing(self): + """ No preprocessing has been done by team DC61 """ + return None + + def get_run_level_analysis(self): + """ No run level analysis has been done by team DC61 """ + return None + + # @staticmethod # Starting python 3.10, staticmethod should be used here + # Otherwise it produces a TypeError: 'staticmethod' object is not callable + def get_subject_information(event_files: list): + """ Create Bunchs for SpecifySPMModel. + + Parameters : + - event_files: list of str, list of events files (one per run) for the subject + + Returns : + - subject_info : list of Bunch for 1st level analysis. + """ + from nipype.interfaces.base import Bunch + + subject_info = [] + + for run_id, event_file in enumerate(event_files): + + onsets = [] + durations = [] + gain_value = [] + loss_value = [] + reaction_time = [] + + # Parse the events file + with open(event_file, 'rt') as file: + next(file) # skip the header + + for line in file: + info = line.strip().split() + + onsets.append(float(info[0])) + durations.append(float(info[1])) + gain_value.append(float(info[2])) + loss_value.append(float(info[3])) + reaction_time.append(float(info[4])) + + # Create a Bunch for the run + subject_info.append( + Bunch( + conditions = [f'gamble_run{run_id + 1}'], + onsets = [onsets], + durations = [durations], + amplitudes = None, + tmod = None, + pmod = [ + Bunch( + name = ['gain_param', 'loss_param', 'rt_param'], + poly = [1, 1, 1], + param = [gain_value, loss_value, reaction_time] + ) + ], + regressor_names = None, + regressors = None + )) + + return subject_info + + # @staticmethod # Starting python 3.10, staticmethod should be used here + # Otherwise it produces a TypeError: 'staticmethod' object is not callable + def get_confounds_file(filepath, subject_id, run_id, working_dir): + """ + Create a new tsv files with only desired confounds per subject per run. + Also computes the first derivative of the motion parameters. + + Parameters : + - filepath : path to the subject confounds file + - subject_id : related subject id + - run_id : related run id + - working_dir: str, name of the directory for intermediate results + + Return : + - confounds_file : path to new file containing only desired confounds + """ + from os import makedirs + from os.path import join + + from pandas import DataFrame, read_csv + from numpy import array, transpose, diff, square, insert + + # Open original confounds file + data_frame = read_csv(filepath, sep = '\t', header=0) + + # Extract confounds we want to use for the model + retained_parameters = DataFrame(transpose(array([ + data_frame['X'], data_frame['Y'], data_frame['Z'], + data_frame['RotX'], data_frame['RotY'], data_frame['RotZ'], + insert(diff(data_frame['X']), 0, 0), + insert(diff(data_frame['Y']), 0, 0), + insert(diff(data_frame['Z']), 0, 0), + insert(diff(data_frame['RotX']), 0, 0), + insert(diff(data_frame['RotY']), 0, 0), + insert(diff(data_frame['RotZ']), 0, 0), + square(data_frame['X']), square(data_frame['Y']), square(data_frame['Z']), + square(data_frame['RotX']), square(data_frame['RotY']), square(data_frame['RotZ']), + insert(square(diff(data_frame['X'])), 0, 0), + insert(square(diff(data_frame['Y'])), 0, 0), + insert(square(diff(data_frame['Z'])), 0, 0), + insert(square(diff(data_frame['RotX'])), 0, 0), + insert(square(diff(data_frame['RotY'])), 0, 0), + insert(square(diff(data_frame['RotZ'])), 0, 0) + ]))) + + # Write confounds to a file + confounds_file = join(working_dir, 'confounds_files', + f'confounds_file_sub-{subject_id}_run-{run_id}.tsv') + + makedirs(join(working_dir, 'confounds_files'), exist_ok = True) + + with open(confounds_file, 'w', encoding = 'utf-8') as writer: + writer.write(retained_parameters.to_csv( + sep = '\t', index = False, header = False, na_rep = '0.0')) + + return confounds_file + + def get_subject_level_analysis(self): + """ + Create the subject level analysis workflow. + + Returns: + - subject_level : nipype.WorkFlow + """ + # Initialize preprocessing workflow to connect nodes along the way + subject_level = Workflow( + base_dir = self.directories.working_dir, name = 'subject_level' + ) + + # Identity interface Node - to iterate over subject_id and run + info_source = Node( + IdentityInterface(fields = ['subject_id']), + name = 'info_source') + info_source.iterables = [('subject_id', self.subject_list)] + + # Select files from derivatives + templates = { + 'func': join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-*_bold_space-MNI152NLin2009cAsym_preproc.nii.gz'), + 'confounds' : join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-*_bold_confounds.tsv'), + 'events': join('sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-*_events.tsv') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory = self.directories.dataset_dir + select_files.inputs.sort_filelist = True + subject_level.connect(info_source, 'subject_id', select_files, 'subject_id') + + # Gunzip - gunzip files because SPM do not use .nii.gz files + gunzip = MapNode(Gunzip(), name = 'gunzip', iterfield=['in_file']) + subject_level.connect(select_files, 'func', gunzip, 'in_file') + + # Smoothing - smooth the func data + smooth = Node(Smooth(), name = 'smooth') + smooth.inputs.fwhm = self.fwhm + smooth.overwrite = False + subject_level.connect(gunzip, 'out_file', smooth, 'in_files') + + # Function node get_subject_info - get subject specific condition information + subject_info = Node(Function( + function = self.get_subject_information, + input_names = ['event_files'], + output_names = ['subject_info'] + ), + name = 'subject_info') + subject_level.connect(select_files, 'events', subject_info, 'event_files') + + # Function node get_confounds_file - Generate confounds files + confounds = MapNode(Function( + function = self.get_confounds_file, + input_names = ['filepath', 'subject_id', 'run_id', 'working_dir'], + output_names = ['confounds_file']), + name = 'confounds', iterfield = ['filepath', 'run_id']) + confounds.inputs.working_dir = self.directories.working_dir + confounds.inputs.run_id = self.run_list + subject_level.connect(info_source, 'subject_id', confounds, 'subject_id') + subject_level.connect(select_files, 'confounds', confounds, 'filepath') + + specify_model = Node(SpecifySPMModel(), name = 'specify_model') + specify_model.inputs.concatenate_runs = False + specify_model.inputs.input_units = 'secs' + specify_model.inputs.output_units = 'secs' + specify_model.inputs.time_repetition = TaskInformation()['RepetitionTime'] + specify_model.inputs.high_pass_filter_cutoff = 128 + specify_model.overwrite = False + subject_level.connect(subject_info, 'subject_info', specify_model, 'subject_info') + subject_level.connect(confounds, 'confounds_file', specify_model, 'realignment_parameters') + subject_level.connect(smooth, 'smoothed_files', specify_model, 'functional_runs') + + model_design = Node(Level1Design(), name = 'model_design') + model_design.inputs.bases = {'hrf': {'derivs': [1, 0]}} # Temporal derivatives + model_design.inputs.timing_units = 'secs' + model_design.inputs.interscan_interval = TaskInformation()['RepetitionTime'] + model_design.overwrite = False + subject_level.connect(specify_model, 'session_info', model_design, 'session_info') + + model_estimate = Node(EstimateModel(), name = 'model_estimate') + model_estimate.inputs.estimation_method = {'Classical': 1} + model_estimate.overwrite = False + subject_level.connect(model_design, 'spm_mat_file', model_estimate, 'spm_mat_file') + + contrast_estimate = Node(EstimateContrast(), name = 'contraste_estimate') + contrast_estimate.inputs.contrasts = self.subject_level_contrasts + contrast_estimate.config = {'execution': {'remove_unnecessary_outputs': False}} + contrast_estimate.overwrite = False + subject_level.connect(model_estimate, 'spm_mat_file', contrast_estimate, 'spm_mat_file') + subject_level.connect(model_estimate, 'beta_images', contrast_estimate, 'beta_images') + subject_level.connect( + model_estimate, 'residual_image', contrast_estimate, 'residual_image') + + # DataSink - store the wanted results in the wanted repository + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + subject_level.connect( + contrast_estimate, 'con_images', data_sink, f'{subject_level.name}.@con_images') + subject_level.connect( + contrast_estimate, 'spm_mat_file', data_sink, f'{subject_level.name}.@spm_mat_file') + + # Remove large files, if requested + if Configuration()['pipelines']['remove_unused_data']: + + # Remove Node - Remove gunzip files once they are no longer needed + remove_gunzip = MapNode( + InterfaceFactory.create('remove_parent_directory'), + name = 'remove_gunzip', + iterfield = ['file_name'] + ) + + # Remove Node - Remove smoothed files once they are no longer needed + remove_smooth = MapNode( + InterfaceFactory.create('remove_parent_directory'), + name = 'remove_smooth', + iterfield = ['file_name'] + ) + + # Add connections + subject_level.connect([ + (smooth, remove_gunzip, [('smoothed_files', '_')]), + (gunzip, remove_gunzip, [('out_file', 'file_name')]), + (data_sink, remove_smooth, [('out_file', '_')]), + (smooth, remove_smooth, [('smoothed_files', 'file_name')]) + ]) + + return subject_level + + def get_subject_level_outputs(self): + """ Return the names of the files the subject level analysis is supposed to generate. """ + + templates = [join( + self.directories.output_dir, + 'subject_level', '_subject_id_{subject_id}', f'con_{contrast_id}.nii')\ + for contrast_id in self.contrast_list] + templates += [join( + self.directories.output_dir, + 'subject_level', '_subject_id_{subject_id}', 'SPM.mat')] + templates += [join( + self.directories.output_dir, + 'subject_level', '_subject_id_{subject_id}', f'spmT_{contrast_id}.nii')\ + for contrast_id in self.contrast_list] + + # Format with subject_ids + return_list = [] + for template in templates: + return_list += [template.format(subject_id = s) for s in self.subject_list] + + return return_list + + def get_group_level_analysis(self): + """ + Return all workflows for the group level analysis. + + Returns; + - a list of nipype.WorkFlow + """ + + return [ + self.get_group_level_analysis_single_group('equalRange'), + self.get_group_level_analysis_single_group('equalIndifference'), + self.get_group_level_analysis_group_comparison() + ] + + def get_group_level_analysis_single_group(self, method): + """ + Return a workflow for the group level analysis in the single group case. + + Parameters: + - method: one of 'equalRange', 'equalIndifference' + + Returns: + - group_level_analysis: nipype.WorkFlow + """ + # Compute the number of participants used to do the analysis + nb_subjects = len(self.subject_list) + + # Infosource - a function free node to iterate over the list of subject names + info_source = Node(IdentityInterface(fields=['contrast_id']), + name = 'info_source') + info_source.iterables = [('contrast_id', self.contrast_list)] + + # Select files from subject level analysis + templates = { + 'contrasts': join(self.directories.output_dir, + 'subject_level', '_subject_id_*', 'con_{contrast_id}.nii'), + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.sort_filelist = True + select_files.inputs.base_directory = self.directories.dataset_dir + + # Datasink - save important files + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + + # Function Node get_group_subjects + # Get subjects in the group and in the subject_list + get_group_subjects = Node(Function( + function = list_intersection, + input_names = ['list_1', 'list_2'], + output_names = ['out_list'] + ), + name = 'get_group_subjects' + ) + get_group_subjects.inputs.list_1 = get_group(method) + get_group_subjects.inputs.list_2 = self.subject_list + + # Create a function to complete the subject ids out from the get_equal_*_subjects nodes + # If not complete, subject id '001' in search patterns + # would match all contrast files with 'con_0001.nii'. + complete_subject_ids = lambda l : [f'_subject_id_{a}' for a in l] + + # Function Node elements_in_string + # Get contrast files for required subjects + # Note : using a MapNode with elements_in_string requires using clean_list to remove + # None values from the out_list + get_contrasts = MapNode(Function( + function = elements_in_string, + input_names = ['input_str', 'elements'], + output_names = ['out_list'] + ), + name = 'get_contrasts', iterfield = 'input_str' + ) + + # One Sample T-Test Design - creates one sample T-Test Design + onesamplettestdes = Node(OneSampleTTestDesign(), name = 'onesampttestdes') + + # EstimateModel - estimate the parameters of the model + # Even for second level it should be 'Classical': 1. + level2estimate = Node(EstimateModel(), name = 'level2estimate') + level2estimate.inputs.estimation_method = {'Classical': 1} + + # EstimateContrast - estimates simple group contrast + level2conestimate = Node(EstimateContrast(), name = 'level2conestimate') + level2conestimate.inputs.group_contrast = True + level2conestimate.inputs.contrasts = [ + ['Group', 'T', ['mean'], [1]], ['Group', 'T', ['mean'], [-1]]] + + # Threshold Node - Create thresholded maps + threshold = MapNode(Threshold(), name = 'threshold', + iterfield = ['stat_image', 'contrast_index']) + threshold.inputs.use_fwe_correction = True + threshold.inputs.height_threshold_type = 'p-value' + threshold.inputs.force_activation = False + threshold.inputs.height_threshold = 0.05 + threshold.inputs.contrast_index = [1, 2] + + # Create the group level workflow + group_level_analysis = Workflow( + base_dir = self.directories.working_dir, + name = f'group_level_analysis_{method}_nsub_{nb_subjects}') + group_level_analysis.connect([ + (info_source, select_files, [('contrast_id', 'contrast_id')]), + (select_files, get_contrasts, [('contrasts', 'input_str')]), + (get_group_subjects, get_contrasts, [ + (('out_list', complete_subject_ids), 'elements') + ]), + (get_contrasts, onesamplettestdes, [ + (('out_list', clean_list), 'in_files') + ]), + #(select_files, onesamplettestdes, [('mask', 'explicit_mask_file')]), + (onesamplettestdes, level2estimate, [('spm_mat_file', 'spm_mat_file')]), + (level2estimate, level2conestimate, [ + ('spm_mat_file', 'spm_mat_file'), + ('beta_images', 'beta_images'), + ('residual_image', 'residual_image') + ]), + (level2conestimate, threshold, [ + ('spm_mat_file', 'spm_mat_file'), + ('spmT_images', 'stat_image') + ]), + (level2estimate, data_sink, [ + ('mask_image', f'{group_level_analysis.name}.@mask')]), + (level2conestimate, data_sink, [ + ('spm_mat_file', f'{group_level_analysis.name}.@spm_mat'), + ('spmT_images', f'{group_level_analysis.name}.@T'), + ('con_images', f'{group_level_analysis.name}.@con')]), + (threshold, data_sink, [ + ('thresholded_map', f'{group_level_analysis.name}.@thresh')]) + ]) + + return group_level_analysis + + def get_group_level_analysis_group_comparison(self): + """ + Return a workflow for the group level analysis in the group comparison case. + + Returns: + - group_level_analysis: nipype.WorkFlow + """ + # Compute the number of participants used to do the analysis + nb_subjects = len(self.subject_list) + + # Infosource - a function free node to iterate over the list of subject names + info_source = Node(IdentityInterface(fields=['contrast_id']), + name = 'info_source') + info_source.iterables = [('contrast_id', self.contrast_list)] + + # Select files from subject level analysis + templates = { + 'contrasts': join(self.directories.output_dir, + 'subject_level', '_subject_id_*', 'con_{contrast_id}.nii'), + #'mask': join('derivatives/fmriprep/gr_mask_tmax.nii') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.sort_filelist = True + select_files.inputs.base_directory = self.directories.dataset_dir + + # Datasink - save important files + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + + # Function Node get_group_subjects + # Get subjects in the group and in the subject_list + get_equal_indifference_subjects = Node(Function( + function = list_intersection, + input_names = ['list_1', 'list_2'], + output_names = ['out_list'] + ), + name = 'get_equal_indifference_subjects' + ) + get_equal_indifference_subjects.inputs.list_1 = get_group('equalIndifference') + get_equal_indifference_subjects.inputs.list_2 = self.subject_list + + # Function Node get_group_subjects + # Get subjects in the group and in the subject_list + get_equal_range_subjects = Node(Function( + function = list_intersection, + input_names = ['list_1', 'list_2'], + output_names = ['out_list'] + ), + name = 'get_equal_range_subjects' + ) + get_equal_range_subjects.inputs.list_1 = get_group('equalRange') + get_equal_range_subjects.inputs.list_2 = self.subject_list + + # Create a function to complete the subject ids out from the get_equal_*_subjects nodes + # If not complete, subject id '001' in search patterns + # would match all contrast files with 'con_0001.nii'. + complete_subject_ids = lambda l : [f'_subject_id_{a}' for a in l] + + # Function Node elements_in_string + # Get contrast files for required subjects + # Note : using a MapNode with elements_in_string requires using clean_list to remove + # None values from the out_list + get_equal_indifference_contrasts = MapNode(Function( + function = elements_in_string, + input_names = ['input_str', 'elements'], + output_names = ['out_list'] + ), + name = 'get_equal_indifference_contrasts', iterfield = 'input_str' + ) + get_equal_range_contrasts = MapNode(Function( + function = elements_in_string, + input_names = ['input_str', 'elements'], + output_names = ['out_list'] + ), + name = 'get_equal_range_contrasts', iterfield = 'input_str' + ) + + # Two Sample T-Test Design + twosampttest = Node(TwoSampleTTestDesign(), name = 'twosampttest') + + # EstimateModel - estimate the parameters of the model + # Even for second level it should be 'Classical': 1. + level2estimate = Node(EstimateModel(), name = 'level2estimate') + level2estimate.inputs.estimation_method = {'Classical': 1} + + # EstimateContrast - estimates simple group contrast + level2conestimate = Node(EstimateContrast(), name = 'level2conestimate') + level2conestimate.inputs.group_contrast = True + level2conestimate.inputs.contrasts = [ + ['Eq range vs Eq indiff in loss', 'T', ['Group_{1}', 'Group_{2}'], [1, -1]] + ] + + # Threshold Node - Create thresholded maps + threshold = Node(Threshold(), name = 'threshold') + threshold.inputs.use_fwe_correction = True + threshold.inputs.height_threshold_type = 'p-value' + threshold.inputs.force_activation = False + threshold.inputs.height_threshold = 0.05 + threshold.inputs.contrast_index = 1 + + # Create the group level workflow + group_level_analysis = Workflow( + base_dir = self.directories.working_dir, + name = f'group_level_analysis_groupComp_nsub_{nb_subjects}') + group_level_analysis.connect([ + (info_source, select_files, [('contrast_id', 'contrast_id')]), + (select_files, get_equal_range_contrasts, [('contrasts', 'input_str')]), + (select_files, get_equal_indifference_contrasts, [('contrasts', 'input_str')]), + (get_equal_range_subjects, get_equal_range_contrasts, [ + (('out_list', complete_subject_ids), 'elements') + ]), + (get_equal_indifference_subjects, get_equal_indifference_contrasts, [ + (('out_list', complete_subject_ids), 'elements') + ]), + (get_equal_range_contrasts, twosampttest, [ + (('out_list', clean_list), 'group1_files') + ]), + (get_equal_indifference_contrasts, twosampttest, [ + (('out_list', clean_list), 'group2_files') + ]), + #(select_files, twosampttest, [('mask', 'explicit_mask_file')]), + (twosampttest, level2estimate, [('spm_mat_file', 'spm_mat_file')]), + (level2estimate, level2conestimate, [ + ('spm_mat_file', 'spm_mat_file'), + ('beta_images', 'beta_images'), + ('residual_image', 'residual_image') + ]), + (level2conestimate, threshold, [ + ('spm_mat_file', 'spm_mat_file'), + ('spmT_images', 'stat_image') + ]), + (level2estimate, data_sink, [ + ('mask_image', f'{group_level_analysis.name}.@mask')]), + (level2conestimate, data_sink, [ + ('spm_mat_file', f'{group_level_analysis.name}.@spm_mat'), + ('spmT_images', f'{group_level_analysis.name}.@T'), + ('con_images', f'{group_level_analysis.name}.@con')]), + (threshold, data_sink, [ + ('thresholded_map', f'{group_level_analysis.name}.@thresh')]) + ]) + + return group_level_analysis + + def get_group_level_outputs(self): + """ Return all names for the files the group level analysis is supposed to generate. """ + + # Handle equalRange and equalIndifference + parameters = { + 'contrast_id': self.contrast_list, + 'method': ['equalRange', 'equalIndifference'], + 'file': [ + 'con_0001.nii', 'con_0002.nii', 'mask.nii', 'SPM.mat', + 'spmT_0001.nii', 'spmT_0002.nii', + join('_threshold0', 'spmT_0001_thr.nii'), join('_threshold1', 'spmT_0002_thr.nii') + ], + 'nb_subjects' : [str(len(self.subject_list))] + } + + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'group_level_analysis_{method}_nsub_{nb_subjects}', + '_contrast_id_{contrast_id}', + '{file}' + ) + return_list = [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + # Handle groupComp + parameters = { + 'contrast_id': self.contrast_list, + 'method': ['groupComp'], + 'file': [ + 'con_0001.nii', 'mask.nii', 'SPM.mat', 'spmT_0001.nii', + join('_threshold0', 'spmT_0001_thr.nii') + ], + 'nb_subjects' : [str(len(self.subject_list))] + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'group_level_analysis_{method}_nsub_{nb_subjects}', + '_contrast_id_{contrast_id}', + '{file}' + ) + return_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + return return_list + + def get_hypotheses_outputs(self): + """ Return all hypotheses output file names. """ + nb_sub = len(self.subject_list) + files = [ + # Hypothesis 1 + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0001', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0001', 'spmT_0001.nii'), + # Hypothesis 2 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0001', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0001', 'spmT_0001.nii'), + # Hypothesis 3 + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0001', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0001', 'spmT_0001.nii'), + # Hypothesis 4 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0001', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0001', 'spmT_0001.nii'), + # Hypothesis 5 + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0002', '_threshold1', 'spmT_0002_thr.nii'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0002', 'spmT_0002.nii'), + # Hypothesis 6 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0002', '_threshold1', 'spmT_0002_thr.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0002', 'spmT_0002.nii'), + # Hypothesis 7 + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0002', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0002', 'spmT_0001.nii'), + # Hypothesis 8 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0002', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0002', 'spmT_0001.nii'), + # Hypothesis 9 + join(f'group_level_analysis_groupComp_nsub_{nb_sub}', + '_contrast_id_0002', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_groupComp_nsub_{nb_sub}', + '_contrast_id_0002', 'spmT_0001.nii') + ] + return [join(self.directories.output_dir, f) for f in files] diff --git a/tests/pipelines/test_team_DC61.py b/tests/pipelines/test_team_DC61.py new file mode 100644 index 00000000..3353b3c3 --- /dev/null +++ b/tests/pipelines/test_team_DC61.py @@ -0,0 +1,128 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Tests of the 'narps_open.pipelines.team_DC61' module. + +Launch this test with PyTest + +Usage: +====== + pytest -q test_team_DC61.py + pytest -q test_team_DC61.py -k +""" +from os.path import join, exists +from filecmp import cmp + +from pytest import helpers, mark +from nipype import Workflow +from nipype.interfaces.base import Bunch + +from narps_open.utils.configuration import Configuration +from narps_open.pipelines.team_DC61 import PipelineTeamDC61 + +class TestPipelinesTeamDC61: + """ A class that contains all the unit tests for the PipelineTeamDC61 class.""" + + @staticmethod + @mark.unit_test + def test_create(): + """ Test the creation of a PipelineTeamDC61 object """ + + pipeline = PipelineTeamDC61() + + # 1 - check the parameters + assert pipeline.fwhm == 5.0 + assert pipeline.team_id == 'DC61' + + # 2 - check workflows + assert pipeline.get_preprocessing() is None + assert pipeline.get_run_level_analysis() is None + assert isinstance(pipeline.get_subject_level_analysis(), Workflow) + group_level = pipeline.get_group_level_analysis() + assert len(group_level) == 3 + for sub_workflow in group_level: + assert isinstance(sub_workflow, Workflow) + + @staticmethod + @mark.unit_test + def test_outputs(): + """ Test the expected outputs of a PipelineTeamDC61 object """ + pipeline = PipelineTeamDC61() + # 1 - 1 subject outputs + pipeline.subject_list = ['001'] + helpers.test_pipeline_outputs(pipeline, [0, 0, 5, 8*2*2 + 5*2, 18]) + + # 2 - 4 subjects outputs + pipeline.subject_list = ['001', '002', '003', '004'] + helpers.test_pipeline_outputs(pipeline, [0, 0, 20, 8*2*2 + 5*2, 18]) + + @staticmethod + @mark.unit_test + def test_subject_information(): + """ Test the get_subject_information method """ + + # Get test files + test_file = join(Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') + info = PipelineTeamDC61.get_subject_information([test_file, test_file]) + + # Compare bunches to expected + bunch = info[0] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['gamble_run1'] + helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 19.535, 27.535, 36.435]]) + helpers.compare_float_2d_arrays(bunch.durations, [[4.0, 4.0, 4.0, 4.0, 4.0]]) + assert bunch.amplitudes is None + assert bunch.tmod is None + assert bunch.pmod[0].name == ['gain_param', 'loss_param', 'rt_param'] + assert bunch.pmod[0].poly == [1, 1, 1] + helpers.compare_float_2d_arrays(bunch.pmod[0].param, [ + [14.0, 34.0, 38.0, 10.0, 16.0], + [6.0, 14.0, 19.0, 15.0, 17.0], + [2.388, 2.289, 0.0, 2.08, 2.288] + ]) + assert bunch.regressor_names is None + assert bunch.regressors is None + + bunch = info[1] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['gamble_run2'] + helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 19.535, 27.535, 36.435]]) + helpers.compare_float_2d_arrays(bunch.durations, [[4.0, 4.0, 4.0, 4.0, 4.0]]) + assert bunch.amplitudes is None + assert bunch.tmod is None + assert bunch.pmod[0].name == ['gain_param', 'loss_param', 'rt_param'] + assert bunch.pmod[0].poly == [1, 1, 1] + helpers.compare_float_2d_arrays(bunch.pmod[0].param, [ + [14.0, 34.0, 38.0, 10.0, 16.0], + [6.0, 14.0, 19.0, 15.0, 17.0], + [2.388, 2.289, 0.0, 2.08, 2.288] + ]) + assert bunch.regressor_names is None + assert bunch.regressors is None + + @staticmethod + @mark.unit_test + def test_confounds_file(temporary_data_dir): + """ Test the get_confounds_file method """ + + confounds_file = join( + Configuration()['directories']['test_data'], 'pipelines', 'confounds.tsv') + reference_file = join( + Configuration()['directories']['test_data'], 'pipelines', 'team_DC61', 'confounds.tsv') + + # Get new confounds file + PipelineTeamDC61.get_confounds_file(confounds_file, 'sid', 'rid', temporary_data_dir) + + # Check confounds file was created + created_confounds_file = join( + temporary_data_dir, 'confounds_files', 'confounds_file_sub-sid_run-rid.tsv') + assert exists(created_confounds_file) + + # Check contents + assert cmp(reference_file, created_confounds_file) + + @staticmethod + @mark.pipeline_test + def test_execution(): + """ Test the execution of a PipelineTeamDC61 and compare results """ + helpers.test_pipeline_evaluation('DC61') diff --git a/tests/test_data/pipelines/team_DC61/confounds.tsv b/tests/test_data/pipelines/team_DC61/confounds.tsv new file mode 100644 index 00000000..1453081d --- /dev/null +++ b/tests/test_data/pipelines/team_DC61/confounds.tsv @@ -0,0 +1,3 @@ +0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +-0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 -0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 9.937996410250001e-05 0.00098247141136 9.0559466761e-12 1.7605839969e-06 1.4760426124900002e-07 2.82878761e-08 9.937996410250001e-05 0.00098247141136 9.0559466761e-12 1.7605839969e-06 1.4760426124900002e-07 2.82878761e-08 +-2.56954e-05 -0.00923735 0.0549667 0.000997278 -0.00019745 -0.000398988 0.009943254600000001 0.022107050000000003 0.05496970931 -0.00032959199999999986 0.000186743 -0.00023079800000000002 6.6025358116e-10 8.53286350225e-05 0.00302133810889 9.94563409284e-07 3.89865025e-08 1.59191424144e-07 9.886831204042119e-05 0.0004887216597025001 0.003021668941625901 1.0863088646399991e-07 3.4872948049e-08 5.326771680400001e-08 From 8c9bb9860f568f6146f625d8f7e4801ba6b0bac9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Wed, 10 Apr 2024 17:42:31 +0200 Subject: [PATCH 02/12] With group level analysis --- narps_open/pipelines/team_DC61.py | 307 +++++++++++++++++------------- tests/pipelines/test_team_DC61.py | 39 +++- 2 files changed, 213 insertions(+), 133 deletions(-) diff --git a/narps_open/pipelines/team_DC61.py b/narps_open/pipelines/team_DC61.py index ba04a626..3c5dae4e 100644 --- a/narps_open/pipelines/team_DC61.py +++ b/narps_open/pipelines/team_DC61.py @@ -71,7 +71,7 @@ def get_subject_information(event_files: list): gain_value = [] loss_value = [] reaction_time = [] - + # Parse the events file with open(event_file, 'rt') as file: next(file) # skip the header @@ -314,6 +314,51 @@ def get_subject_level_outputs(self): return return_list + def get_group_level_contrasts(subject_level_contrast: str): + """ Return a contrast list for EstimateContrast + + Parameters : + - subject_level_contrast: str, id of the subject level contrast : + either 'effect_of_gain' + or 'effect_of_loss' + + Returns : + - list of contrasts for the EstimateContrast interface. + """ + if subject_level_contrast == 'effect_of_gain': + return [ + ['gain_param_range', 'T', ['equalIndifference', 'equalRange'], [0, 1]], + ['gain_param_indiff', 'T', ['equalIndifference', 'equalRange'], [1, 0]] + ] + + if subject_level_contrast == 'effect_of_loss': + return [ + ['loss_param_range', 'F', ['equalIndifference', 'equalRange'], [0, 1]], + ['loss_param_indiff', 'F', ['equalIndifference', 'equalRange'], [1, 0]] + ] + + return [] + + def get_group_covariates(subjects: list): + """ Return a covariates list for OneSampleTTestDesign + + Parameters : + - subjects: list, list of the subjects involved in the design + + Returns : + - list of covariates for the OneSampleTTestDesign interface. + """ + from narps_open.data.participants import get_group + + return [ + dict( + vector = [1 if s in get_group('equalRange') else 0 for s in subjects], + name = 'equalRange'), + dict( + vector = [1 if s in get_group('equalIndifference') else 0 for s in subjects], + name = 'equalIndifference') + ] + def get_group_level_analysis(self): """ Return all workflows for the group level analysis. @@ -321,88 +366,106 @@ def get_group_level_analysis(self): Returns; - a list of nipype.WorkFlow """ - return [ - self.get_group_level_analysis_single_group('equalRange'), - self.get_group_level_analysis_single_group('equalIndifference'), + self.get_group_level_analysis_single_group(), self.get_group_level_analysis_group_comparison() - ] + ] - def get_group_level_analysis_single_group(self, method): + def get_group_level_analysis_single_group(self): """ - Return a workflow for the group level analysis in the single group case. - - Parameters: - - method: one of 'equalRange', 'equalIndifference' + Return a workflow for the group level analysis. Returns: - - group_level_analysis: nipype.WorkFlow + - group_level: nipype.WorkFlow """ # Compute the number of participants used to do the analysis nb_subjects = len(self.subject_list) - # Infosource - a function free node to iterate over the list of subject names - info_source = Node(IdentityInterface(fields=['contrast_id']), + # Create the group level workflow + group_level = Workflow( + base_dir = self.directories.working_dir, + name = f'group_level_analysis_nsub_{nb_subjects}') + + # IDENTITY INTERFACE - Iterate over the list of subject-level contrasts + info_source = Node(IdentityInterface(fields=['contrast_id', 'contrast_name']), name = 'info_source') - info_source.iterables = [('contrast_id', self.contrast_list)] + info_source.iterables = [ + ('contrast_id', self.contrast_list), + ('contrast_name', [c[0] for c in self.subject_level_contrasts])] + info_source.synchronize = True - # Select files from subject level analysis + # SELECT FILES - Get files from subject-level analysis templates = { 'contrasts': join(self.directories.output_dir, - 'subject_level', '_subject_id_*', 'con_{contrast_id}.nii'), + 'subject_level', '_subject_id_*', 'con_{contrast_id}.nii') } select_files = Node(SelectFiles(templates), name = 'select_files') select_files.inputs.sort_filelist = True select_files.inputs.base_directory = self.directories.dataset_dir - # Datasink - save important files - data_sink = Node(DataSink(), name = 'data_sink') - data_sink.inputs.base_directory = self.directories.output_dir - - # Function Node get_group_subjects - # Get subjects in the group and in the subject_list - get_group_subjects = Node(Function( - function = list_intersection, - input_names = ['list_1', 'list_2'], - output_names = ['out_list'] - ), - name = 'get_group_subjects' - ) - get_group_subjects.inputs.list_1 = get_group(method) - get_group_subjects.inputs.list_2 = self.subject_list - # Create a function to complete the subject ids out from the get_equal_*_subjects nodes # If not complete, subject id '001' in search patterns # would match all contrast files with 'con_0001.nii'. complete_subject_ids = lambda l : [f'_subject_id_{a}' for a in l] + # GET CONTRAST FILES # Function Node elements_in_string # Get contrast files for required subjects # Note : using a MapNode with elements_in_string requires using clean_list to remove # None values from the out_list - get_contrasts = MapNode(Function( + get_contrast_files = MapNode(Function( function = elements_in_string, input_names = ['input_str', 'elements'], output_names = ['out_list'] ), - name = 'get_contrasts', iterfield = 'input_str' + name = 'get_contrast_files', iterfield = 'input_str' + ) + get_contrast_files.inputs.elements = complete_subject_ids(self.subject_list) + group_level.connect(select_files, 'contrasts', get_contrast_files, 'input_str') + + # GET COVARIATES + # Function Node get_group_covariates + # Get groups as covariates input for OneSampleTTestDesign + get_covariates = Node(Function( + function = self.get_group_covariates, + input_names = ['subjects'], + output_names = ['covariates'] + ), + name = 'get_covariates' ) + get_covariates.inputs.subjects = self.subject_list - # One Sample T-Test Design - creates one sample T-Test Design - onesamplettestdes = Node(OneSampleTTestDesign(), name = 'onesampttestdes') + # ONE SAMPLE T-TEST DESIGN - Create a t test design for hypothesis testing inside a group + one_sample_t_test = Node(OneSampleTTestDesign(), name = 'one_sample_t_test') + group_level.connect(get_covariates, 'covariates', one_sample_t_test, 'covariates') + group_level.connect( + get_contrast_files, ('out_list', clean_list), one_sample_t_test, 'in_files') - # EstimateModel - estimate the parameters of the model + # ESTIMATE MODEL - estimate the parameters of the model # Even for second level it should be 'Classical': 1. - level2estimate = Node(EstimateModel(), name = 'level2estimate') - level2estimate.inputs.estimation_method = {'Classical': 1} + model_estimate = Node(EstimateModel(), name = 'model_estimate') + model_estimate.inputs.estimation_method = {'Classical': 1} + group_level.connect(one_sample_t_test, 'spm_mat_file', model_estimate, 'spm_mat_file') + + # Function Node get_group_level_contrasts + # Get a contrast list for EstimateContrast + get_contrasts = Node(Function( + function = self.get_group_level_contrasts, + input_names = ['subject_level_contrast'], + output_names = ['contrasts'] + ), + name = 'get_contrasts' + ) - # EstimateContrast - estimates simple group contrast - level2conestimate = Node(EstimateContrast(), name = 'level2conestimate') - level2conestimate.inputs.group_contrast = True - level2conestimate.inputs.contrasts = [ - ['Group', 'T', ['mean'], [1]], ['Group', 'T', ['mean'], [-1]]] + # ESTIMATE CONTRASTS - estimates simple group contrast + contrast_estimate = Node(EstimateContrast(), name = 'contrast_estimate') + contrast_estimate.inputs.group_contrast = True + group_level.connect(get_contrasts, 'contrasts', contrast_estimate, 'contrasts') + group_level.connect(model_estimate, 'spm_mat_file', contrast_estimate, 'spm_mat_file') + group_level.connect(model_estimate, 'beta_images', contrast_estimate, 'beta_images') + group_level.connect(model_estimate, 'residual_image', contrast_estimate, 'residual_image') - # Threshold Node - Create thresholded maps + # THRESHOLD - Create thresholded maps threshold = MapNode(Threshold(), name = 'threshold', iterfield = ['stat_image', 'contrast_index']) threshold.inputs.use_fwe_correction = True @@ -410,42 +473,24 @@ def get_group_level_analysis_single_group(self, method): threshold.inputs.force_activation = False threshold.inputs.height_threshold = 0.05 threshold.inputs.contrast_index = [1, 2] + group_level.connect(contrast_estimate, 'spm_mat_file', threshold, 'spm_mat_file') + group_level.connect(contrast_estimate, 'spmT_images', threshold, 'stat_image') - # Create the group level workflow - group_level_analysis = Workflow( - base_dir = self.directories.working_dir, - name = f'group_level_analysis_{method}_nsub_{nb_subjects}') - group_level_analysis.connect([ - (info_source, select_files, [('contrast_id', 'contrast_id')]), - (select_files, get_contrasts, [('contrasts', 'input_str')]), - (get_group_subjects, get_contrasts, [ - (('out_list', complete_subject_ids), 'elements') - ]), - (get_contrasts, onesamplettestdes, [ - (('out_list', clean_list), 'in_files') - ]), - #(select_files, onesamplettestdes, [('mask', 'explicit_mask_file')]), - (onesamplettestdes, level2estimate, [('spm_mat_file', 'spm_mat_file')]), - (level2estimate, level2conestimate, [ - ('spm_mat_file', 'spm_mat_file'), - ('beta_images', 'beta_images'), - ('residual_image', 'residual_image') - ]), - (level2conestimate, threshold, [ - ('spm_mat_file', 'spm_mat_file'), - ('spmT_images', 'stat_image') - ]), - (level2estimate, data_sink, [ - ('mask_image', f'{group_level_analysis.name}.@mask')]), - (level2conestimate, data_sink, [ - ('spm_mat_file', f'{group_level_analysis.name}.@spm_mat'), - ('spmT_images', f'{group_level_analysis.name}.@T'), - ('con_images', f'{group_level_analysis.name}.@con')]), + # Datasink - save important files + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + group_level.connect([ + (model_estimate, data_sink, [ + ('mask_image', f'{group_level.name}.@mask')]), + (contrast_estimate, data_sink, [ + ('spm_mat_file', f'{group_level.name}.@spm_mat'), + ('spmT_images', f'{group_level.name}.@T'), + ('con_images', f'{group_level.name}.@con')]), (threshold, data_sink, [ - ('thresholded_map', f'{group_level_analysis.name}.@thresh')]) + ('thresholded_map', f'{group_level.name}.@thresh')]) ]) - return group_level_analysis + return group_level def get_group_level_analysis_group_comparison(self): """ @@ -457,25 +502,21 @@ def get_group_level_analysis_group_comparison(self): # Compute the number of participants used to do the analysis nb_subjects = len(self.subject_list) - # Infosource - a function free node to iterate over the list of subject names - info_source = Node(IdentityInterface(fields=['contrast_id']), - name = 'info_source') - info_source.iterables = [('contrast_id', self.contrast_list)] + # Create the group level workflow + group_level_analysis = Workflow( + base_dir = self.directories.working_dir, + name = f'group_level_analysis_groupComp_nsub_{nb_subjects}') - # Select files from subject level analysis + # SELEC FILES - Get files from subject level analysis (effect_of_loss contrast only) templates = { 'contrasts': join(self.directories.output_dir, - 'subject_level', '_subject_id_*', 'con_{contrast_id}.nii'), - #'mask': join('derivatives/fmriprep/gr_mask_tmax.nii') + 'subject_level', '_subject_id_*', 'con_0002.nii') } select_files = Node(SelectFiles(templates), name = 'select_files') select_files.inputs.sort_filelist = True select_files.inputs.base_directory = self.directories.dataset_dir - # Datasink - save important files - data_sink = Node(DataSink(), name = 'data_sink') - data_sink.inputs.base_directory = self.directories.output_dir - + # GET SUBJECT LIST IN EACH GROUP # Function Node get_group_subjects # Get subjects in the group and in the subject_list get_equal_indifference_subjects = Node(Function( @@ -505,6 +546,7 @@ def get_group_level_analysis_group_comparison(self): # would match all contrast files with 'con_0001.nii'. complete_subject_ids = lambda l : [f'_subject_id_{a}' for a in l] + # GET CONTRAST FILES # Function Node elements_in_string # Get contrast files for required subjects # Note : using a MapNode with elements_in_string requires using clean_list to remove @@ -516,6 +558,12 @@ def get_group_level_analysis_group_comparison(self): ), name = 'get_equal_indifference_contrasts', iterfield = 'input_str' ) + group_level_analysis.connect( + select_files, 'contrasts', get_equal_indifference_contrasts, 'input_str') + group_level_analysis.connect( + get_equal_indifference_subjects, ('out_list', complete_subject_ids), + get_equal_indifference_contrasts, 'elements') + get_equal_range_contrasts = MapNode(Function( function = elements_in_string, input_names = ['input_str', 'elements'], @@ -523,64 +571,61 @@ def get_group_level_analysis_group_comparison(self): ), name = 'get_equal_range_contrasts', iterfield = 'input_str' ) - - # Two Sample T-Test Design - twosampttest = Node(TwoSampleTTestDesign(), name = 'twosampttest') - - # EstimateModel - estimate the parameters of the model + group_level_analysis.connect( + select_files, 'contrasts', get_equal_range_contrasts, 'input_str') + group_level_analysis.connect( + get_equal_range_subjects, ('out_list', complete_subject_ids), + get_equal_range_contrasts, 'elements') + + # TWO SAMPLE T-TEST DESIGN - Create a t test design for group comparison + two_sample_t_test = Node(TwoSampleTTestDesign(), name = 'two_sample_t_test') + group_level_analysis.connect( + get_equal_range_contrasts, ('out_list', clean_list), + two_sample_t_test, 'group1_files') + group_level_analysis.connect( + get_equal_indifference_contrasts, ('out_list', clean_list), + two_sample_t_test, 'group2_files') + + # ESTIMATE MODEL - Estimate the parameters of the model # Even for second level it should be 'Classical': 1. - level2estimate = Node(EstimateModel(), name = 'level2estimate') - level2estimate.inputs.estimation_method = {'Classical': 1} + model_estimate = Node(EstimateModel(), name = 'model_estimate') + model_estimate.inputs.estimation_method = {'Classical': 1} + group_level_analysis.connect( + two_sample_t_test, 'spm_mat_file', model_estimate, 'spm_mat_file') - # EstimateContrast - estimates simple group contrast - level2conestimate = Node(EstimateContrast(), name = 'level2conestimate') - level2conestimate.inputs.group_contrast = True - level2conestimate.inputs.contrasts = [ + # ESTIMATE CONTRASTS - Estimates contrasts + contrast_estimate = Node(EstimateContrast(), name = 'contrast_estimate') + contrast_estimate.inputs.group_contrast = True + contrast_estimate.inputs.contrasts = [ ['Eq range vs Eq indiff in loss', 'T', ['Group_{1}', 'Group_{2}'], [1, -1]] ] + group_level_analysis.connect([ + (model_estimate, contrast_estimate, [ + ('spm_mat_file', 'spm_mat_file'), + ('beta_images', 'beta_images'), + ('residual_image', 'residual_image') + ])]) - # Threshold Node - Create thresholded maps + # THRESHOLD - Create thresholded maps threshold = Node(Threshold(), name = 'threshold') threshold.inputs.use_fwe_correction = True threshold.inputs.height_threshold_type = 'p-value' threshold.inputs.force_activation = False threshold.inputs.height_threshold = 0.05 threshold.inputs.contrast_index = 1 - - # Create the group level workflow - group_level_analysis = Workflow( - base_dir = self.directories.working_dir, - name = f'group_level_analysis_groupComp_nsub_{nb_subjects}') group_level_analysis.connect([ - (info_source, select_files, [('contrast_id', 'contrast_id')]), - (select_files, get_equal_range_contrasts, [('contrasts', 'input_str')]), - (select_files, get_equal_indifference_contrasts, [('contrasts', 'input_str')]), - (get_equal_range_subjects, get_equal_range_contrasts, [ - (('out_list', complete_subject_ids), 'elements') - ]), - (get_equal_indifference_subjects, get_equal_indifference_contrasts, [ - (('out_list', complete_subject_ids), 'elements') - ]), - (get_equal_range_contrasts, twosampttest, [ - (('out_list', clean_list), 'group1_files') - ]), - (get_equal_indifference_contrasts, twosampttest, [ - (('out_list', clean_list), 'group2_files') - ]), - #(select_files, twosampttest, [('mask', 'explicit_mask_file')]), - (twosampttest, level2estimate, [('spm_mat_file', 'spm_mat_file')]), - (level2estimate, level2conestimate, [ - ('spm_mat_file', 'spm_mat_file'), - ('beta_images', 'beta_images'), - ('residual_image', 'residual_image') - ]), - (level2conestimate, threshold, [ + (contrast_estimate, threshold, [ ('spm_mat_file', 'spm_mat_file'), ('spmT_images', 'stat_image') - ]), - (level2estimate, data_sink, [ + ])]) + + # DATA SINK - save important files + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + group_level_analysis.connect([ + (model_estimate, data_sink, [ ('mask_image', f'{group_level_analysis.name}.@mask')]), - (level2conestimate, data_sink, [ + (contrast_estimate, data_sink, [ ('spm_mat_file', f'{group_level_analysis.name}.@spm_mat'), ('spmT_images', f'{group_level_analysis.name}.@T'), ('con_images', f'{group_level_analysis.name}.@con')]), diff --git a/tests/pipelines/test_team_DC61.py b/tests/pipelines/test_team_DC61.py index 3353b3c3..d25a81e3 100644 --- a/tests/pipelines/test_team_DC61.py +++ b/tests/pipelines/test_team_DC61.py @@ -13,13 +13,28 @@ from os.path import join, exists from filecmp import cmp -from pytest import helpers, mark +from pytest import helpers, mark, fixture from nipype import Workflow from nipype.interfaces.base import Bunch from narps_open.utils.configuration import Configuration from narps_open.pipelines.team_DC61 import PipelineTeamDC61 +@fixture +def mock_participants_data(mocker): + """ A fixture to provide mocked data from the test_data directory """ + + mocker.patch( + 'narps_open.data.participants.Configuration', + return_value = { + 'directories': { + 'dataset': join( + Configuration()['directories']['test_data'], + 'data', 'participants') + } + } + ) + class TestPipelinesTeamDC61: """ A class that contains all the unit tests for the PipelineTeamDC61 class.""" @@ -39,7 +54,7 @@ def test_create(): assert pipeline.get_run_level_analysis() is None assert isinstance(pipeline.get_subject_level_analysis(), Workflow) group_level = pipeline.get_group_level_analysis() - assert len(group_level) == 3 + assert len(group_level) == 2 for sub_workflow in group_level: assert isinstance(sub_workflow, Workflow) @@ -121,6 +136,26 @@ def test_confounds_file(temporary_data_dir): # Check contents assert cmp(reference_file, created_confounds_file) + @staticmethod + @mark.unit_test + def test_group_level_contrasts(): + """ Test the get_group_contrasts method """ + + @staticmethod + @mark.unit_test + def test_group_covariates(mock_participants_data): + """ Test the get_group_covariates method """ + subjects = ['001', '002', '003', '004'] + assert PipelineTeamDC61.get_group_covariates(subjects) == [ + dict(vector = [0, 1, 0, 1], name = 'equalRange'), + dict(vector = [1, 0, 1, 0], name = 'equalIndifference') + ] + subjects = ['001', '003', '002', '004'] + assert PipelineTeamDC61.get_group_covariates(subjects) == [ + dict(vector = [0, 0, 1, 1], name = 'equalRange'), + dict(vector = [1, 1, 0, 0], name = 'equalIndifference') + ] + @staticmethod @mark.pipeline_test def test_execution(): From 5100256adef87017610c97b2990062399e43729a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Thu, 11 Apr 2024 16:48:41 +0200 Subject: [PATCH 03/12] Starting B5I6 repro [skip ci] --- narps_open/pipelines/__init__.py | 2 +- narps_open/pipelines/team_B5I6.py | 835 ++++++++++++++++++ tests/pipelines/test_team_B5I6.py | 187 ++++ .../team_B5I6/confounds_with_outliers.tsv | 20 + .../team_B5I6/out_confounds_no_outliers.tsv | 3 + .../team_B5I6/out_confounds_outliers.tsv | 19 + 6 files changed, 1065 insertions(+), 1 deletion(-) create mode 100644 narps_open/pipelines/team_B5I6.py create mode 100644 tests/pipelines/test_team_B5I6.py create mode 100644 tests/test_data/pipelines/team_B5I6/confounds_with_outliers.tsv create mode 100644 tests/test_data/pipelines/team_B5I6/out_confounds_no_outliers.tsv create mode 100644 tests/test_data/pipelines/team_B5I6/out_confounds_outliers.tsv diff --git a/narps_open/pipelines/__init__.py b/narps_open/pipelines/__init__.py index dd6f5a38..e02d2f22 100644 --- a/narps_open/pipelines/__init__.py +++ b/narps_open/pipelines/__init__.py @@ -41,7 +41,7 @@ '9U7M': None, 'AO86': None, 'B23O': None, - 'B5I6': None, + 'B5I6': 'PipelineTeamB5I6', 'C22U': None, 'C88N': 'PipelineTeamC88N', 'DC61': 'PipelineTeamDC61', diff --git a/narps_open/pipelines/team_B5I6.py b/narps_open/pipelines/team_B5I6.py new file mode 100644 index 00000000..ad201a9c --- /dev/null +++ b/narps_open/pipelines/team_B5I6.py @@ -0,0 +1,835 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Write the work of NARPS team B5I6 using Nipype """ + +from os.path import join +from itertools import product + +from nipype import Workflow, Node, MapNode +from nipype.interfaces.utility import IdentityInterface, Function, Split +from nipype.interfaces.io import SelectFiles, DataSink +from nipype.interfaces.fsl import ( + IsotropicSmooth, Level1Design, FEATModel, + L2Model, Merge, FLAMEO, FILMGLS, MultipleRegressDesign, + Cluster, SmoothEstimate, FSLCommand, Randomise + ) +from nipype.algorithms.modelgen import SpecifyModel +from nipype.interfaces.fsl.maths import MultiImageMaths + +from narps_open.utils.configuration import Configuration +from narps_open.pipelines import Pipeline +from narps_open.data.task import TaskInformation +from narps_open.data.participants import get_group +from narps_open.core.common import list_intersection, elements_in_string, clean_list +from narps_open.core.interfaces import InterfaceFactory + +# Setup FSL +FSLCommand.set_default_output_type('NIFTI_GZ') + +class PipelineTeamB5I6(Pipeline): + """ A class that defines the pipeline of team B5I6 """ + + def __init__(self): + super().__init__() + self.fwhm = 5.0 + self.team_id = 'B5I6' + self.contrast_list = ['1', '2', '3', '4', '5', '6'] + self.run_level_contrasts = [ + ('positive_main_effect_of_decision', 'T', ['trial', 'gain', 'loss'], [1, 0, 0]), + ('negative_main_effect_of_decision', 'T', ['trial', 'gain', 'loss'], [-1, 0, 0]), + ('positive_parametric_effect_of_gain', 'T', ['trial', 'gain', 'loss'], [0, 1, 0]), + ('negative_parametric_effect_of_gain', 'T', ['trial', 'gain', 'loss'], [0, -1, 0]), + ('positive_parametric_effect_of_loss', 'T', ['trial', 'gain', 'loss'], [0, 0, 1]), + ('negative_parametric_effect_of_loss', 'T', ['trial', 'gain', 'loss'], [0, 0, -1]) + ] + + def get_preprocessing(self): + """ No preprocessing has been done by team B5I6 """ + return None + + def get_subject_information(event_file): + """ + Create Bunchs for specifyModel. + + Parameters : + - event_file : str, file corresponding to the run and the subject to analyze + + Returns : + - subject_info : list of Bunch for 1st level analysis. + """ + from numpy import mean + from nipype.interfaces.base import Bunch + + onsets_trial = [] + onsets_missed = [] + durations_trial = [] + durations_missed = [] + weights_gain = [] + weights_loss = [] + + with open(event_file, 'rt') as file: + next(file) # skip the header + + for line in file: + info = line.strip().split() + + if 'NoResp' in info[5]: + onsets_missed.append(float(info[0])) + durations_missed.append(float(info[1])) + else: + onsets_trial.append(float(info[0])) + durations_trial.append(float(info[1])) + weights_gain.append(float(info[2])) + weights_loss.append(float(info[3])) + + # Mean center magnitudes + weights_gain = weights_gain - mean(weights_gain) + weights_loss = weights_loss - mean(weights_loss) + + # Add missed trials if any + conditions = ['trial'] + onsets = [onsets_trial] + durations = [durations_trial] + if onsets_missed: + conditions.append('missed') + onsets.append(onsets_missed) + durations.append(durations_missed) + + return [ + Bunch( + conditions = conditions, + onsets = onsets, + durations = durations, + amplitudes = None, + tmod = None, + pmod = [ + Bunch( + name = ['gain', 'loss'], + poly = [1, 1], + param = [weights_gain, weights_loss] + ) + ], + regressor_names = None, + regressors = None + ) + ] + + def get_confounds_file(filepath, subject_id, run_id): + """ + Create a tsv file with only desired confounds per subject per run. + + Parameters : + - filepath : path to the subject confounds file (i.e. one per run) + - subject_id : subject for whom the 1st level analysis is made + - run_id: run for which the 1st level analysis is made + + Return : + - confounds_file : paths to new files containing only desired confounds. + """ + from os.path import abspath + + from pandas import read_csv, DataFrame + from numpy import array, transpose, insert, diff, square, nanpercentile, c_ + + data_frame = read_csv(filepath, sep = '\t', header=0) + motion_regressors = array([ + # Motion regressors + data_frame['X'], data_frame['Y'], data_frame['Z'], + data_frame['RotX'], data_frame['RotY'], data_frame['RotZ'], + # Derivatives + insert(diff(data_frame['X']), 0, 0), + insert(diff(data_frame['Y']), 0, 0), + insert(diff(data_frame['Z']), 0, 0), + insert(diff(data_frame['RotX']), 0, 0), + insert(diff(data_frame['RotY']), 0, 0), + insert(diff(data_frame['RotZ']), 0, 0), + # Squared motion regressors + square(data_frame['X']), square(data_frame['Y']), square(data_frame['Z']), + square(data_frame['RotX']), square(data_frame['RotY']), square(data_frame['RotZ']), + # Squared derivatives + insert(square(diff(data_frame['X'])), 0, 0), + insert(square(diff(data_frame['Y'])), 0, 0), + insert(square(diff(data_frame['Z'])), 0, 0), + insert(square(diff(data_frame['RotX'])), 0, 0), + insert(square(diff(data_frame['RotY'])), 0, 0), + insert(square(diff(data_frame['RotZ'])), 0, 0), + ]) + + # Compute outliers of `non-stdDVARS` + dvars_outliers = array(data_frame['non-stdDVARS']) + print(dvars_outliers) + + percentile_75 = nanpercentile(dvars_outliers, 75) + interquartile_range = percentile_75 - nanpercentile(dvars_outliers, 25) + dvars_outliers = (dvars_outliers > (percentile_75 + 1.5 * interquartile_range)).astype(int) + + print(dvars_outliers) + print(percentile_75 + 1.5 * interquartile_range) + + # Compute outliers of `FramewiseDisplacement` + fd_outliers = array(data_frame['FramewiseDisplacement']) + print(fd_outliers) + + percentile_75 = nanpercentile(fd_outliers, 75) + interquartile_range = percentile_75 - nanpercentile(fd_outliers, 25) + fd_outliers = (fd_outliers > (percentile_75 + 1.5 * interquartile_range)).astype(int) + + print(fd_outliers) + print(percentile_75 + 1.5 * interquartile_range) + + # Build an outlier regressor + outlier_regressor = fd_outliers | dvars_outliers + + # Add the regressor if any outliers in the run + if outlier_regressor.any(): + retained_confounds = DataFrame(c_[transpose(motion_regressors), outlier_regressor]) + else : + retained_confounds = DataFrame(transpose(motion_regressors)) + + # Write confounds to a file + confounds_file = abspath(f'confounds_file_sub-{subject_id}_run-{run_id}.tsv') + with open(confounds_file, 'w', encoding = 'utf-8') as writer: + writer.write(retained_confounds.to_csv( + sep = '\t', index = False, header = False, na_rep = '0.0')) + + return confounds_file + + def get_run_level_analysis(self): + """ + Create the run level analysis workflow. + + Returns: + - run_level : nipype.WorkFlow + """ + # Create run level analysis workflow and connect its nodes + run_level = Workflow( + base_dir = self.directories.working_dir, + name = 'run_level_analysis' + ) + + # IdentityInterface Node - Iterate on subject and runs + information_source = Node(IdentityInterface( + fields = ['subject_id', 'run_id']), + name = 'information_source') + information_source.iterables = [ + ('subject_id', self.subject_list), + ('run_id', self.run_list) + ] + + # SelectFiles - Get necessary files + templates = { + 'func' : join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-{run_id}_bold_space-MNI152NLin2009cAsym_preproc.nii.gz'), + 'events' : join('sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-{run_id}_events.tsv'), + 'confounds' : join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-{run_id}_bold_confounds.tsv') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory = self.directories.dataset_dir + run_level.connect(information_source, 'subject_id', select_files, 'subject_id') + run_level.connect(information_source, 'run_id', select_files, 'run_id') + + # IsotropicSmooth Node - Smoothing data + smoothing_func = Node(IsotropicSmooth(), name = 'smoothing_func') + smoothing_func.inputs.fwhm = self.fwhm + run_level.connect(select_files, 'func', smoothing_func, 'in_file') + + # Get Subject Info - get subject specific condition information + subject_information = Node(Function( + function = self.get_subject_information, + input_names = ['event_file'], + output_names = ['subject_info'] + ), name = 'subject_information') + run_level.connect(select_files, 'events', subject_information, 'event_file') + + # Function Node get_confounds_file - Get files with movement confounds + confounds = Node(Function( + function = self.get_confounds_file, + input_names = ['filepath', 'subject_id', 'run_id'], + output_names = ['confounds_file']), + name = 'confounds') + run_level.connect(information_source, 'subject_id', confounds, 'subject_id') + run_level.connect(information_source, 'run_id', confounds, 'run_id') + run_level.connect(select_files, 'confounds', confounds, 'filepath') + + # SpecifyModel Node - Generate run level model + specify_model = Node(SpecifyModel(), name = 'specify_model') + specify_model.inputs.high_pass_filter_cutoff = 100 + specify_model.inputs.input_units = 'secs' + specify_model.inputs.time_repetition = TaskInformation()['RepetitionTime'] + run_level.connect(confounds, 'confounds_file', specify_model, 'realignment_parameters') + run_level.connect(smoothing_func, 'out_file', specify_model, 'functional_runs') + run_level.connect(subject_information, 'subject_info', specify_model, 'subject_info') + + # Level1Design Node - Generate files for run level computation + model_design = Node(Level1Design(), name = 'model_design') + model_design.inputs.bases = {'dgamma' : {'derivs' : False }} + model_design.inputs.interscan_interval = TaskInformation()['RepetitionTime'] + model_design.inputs.model_serial_correlations = True + model_design.inputs.contrasts = self.run_level_contrasts + run_level.connect(specify_model, 'session_info', model_design, 'session_info') + + # FEATModel Node - Generate run level model + model_generation = Node(FEATModel(), name = 'model_generation') + run_level.connect(model_design, 'ev_files', model_generation, 'ev_files') + run_level.connect(model_design, 'fsf_files', model_generation, 'fsf_file') + + # FILMGLS Node - Estimate first level model + model_estimate = Node(FILMGLS(), name='model_estimate') + run_level.connect(smoothing_func, 'out_file', model_estimate, 'in_file') + run_level.connect(model_generation, 'con_file', model_estimate, 'tcon_file') + run_level.connect(model_generation, 'design_file', model_estimate, 'design_file') + + # DataSink Node - store the wanted results in the wanted directory + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + run_level.connect(model_estimate, 'results_dir', data_sink, 'run_level_analysis.@results') + run_level.connect(model_generation, 'design_file', data_sink, 'run_level_analysis.@design_file') + run_level.connect(model_generation, 'design_image', data_sink, 'run_level_analysis.@design_img') + + # Remove large files, if requested + if Configuration()['pipelines']['remove_unused_data']: + remove_smooth = Node( + InterfaceFactory.create('remove_parent_directory'), + name = 'remove_smooth') + run_level.connect(data_sink, 'out_file', remove_smooth, '_') + run_level.connect(smoothing_func, 'out_file', remove_smooth, 'file_name') + + return run_level + + def get_run_level_outputs(self): + """ Return the names of the files the run level analysis is supposed to generate. """ + + parameters = { + 'run_id' : self.run_list, + 'subject_id' : self.subject_list, + 'contrast_id' : self.contrast_list, + } + parameter_sets = product(*parameters.values()) + output_dir = join(self.directories.output_dir, + 'run_level_analysis', '_run_id_{run_id}_subject_id_{subject_id}') + templates = [ + join(output_dir, 'results', 'cope{contrast_id}.nii.gz'), + join(output_dir, 'results', 'tstat{contrast_id}.nii.gz'), + join(output_dir, 'results', 'varcope{contrast_id}.nii.gz'), + join(output_dir, 'results', 'zstat{contrast_id}.nii.gz') + ] + return [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets for template in templates] + + def get_subject_level_analysis(self): + """ + Create the subject level analysis workflow. + + Returns: + - subject_level_analysis : nipype.WorkFlow + """ + # Second level (single-subject, mean of all four scans) analysis workflow. + subject_level = Workflow( + base_dir = self.directories.working_dir, + name = 'subject_level_analysis') + + # Infosource Node - To iterate on subject and runs + information_source = Node(IdentityInterface( + fields = ['subject_id', 'contrast_id']), + name = 'information_source') + information_source.iterables = [ + ('subject_id', self.subject_list), + ('contrast_id', self.contrast_list) + ] + + # SelectFiles node - to select necessary files + templates = { + 'cope' : join(self.directories.output_dir, + 'run_level_analysis', '_run_id_*_subject_id_{subject_id}', 'results', + 'cope{contrast_id}.nii.gz'), + 'varcope' : join(self.directories.output_dir, + 'run_level_analysis', '_run_id_*_subject_id_{subject_id}', 'results', + 'varcope{contrast_id}.nii.gz'), + 'masks' : join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-{run_id}_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory= self.directories.results_dir + subject_level.connect(information_source, 'subject_id', select_files, 'subject_id') + subject_level.connect(information_source, 'contrast_id', select_files, 'contrast_id') + + # Merge Node - Merge copes files for each subject + merge_copes = Node(Merge(), name = 'merge_copes') + merge_copes.inputs.dimension = 't' + subject_level.connect(select_files, 'cope', merge_copes, 'in_files') + + # Merge Node - Merge varcopes files for each subject + merge_varcopes = Node(Merge(), name = 'merge_varcopes') + merge_varcopes.inputs.dimension = 't' + subject_level.connect(select_files, 'varcope', merge_varcopes, 'in_files') + + # Split Node - Split mask list to serve them as inputs of the MultiImageMaths node. + split_masks = Node(Split(), name = 'split_masks') + split_masks.inputs.splits = [1, len(self.run_list) - 1] + split_masks.inputs.squeeze = True # Unfold one-element splits removing the list + subject_level.connect(select_files, 'masks', split_masks, 'inlist') + + # MultiImageMaths Node - Create a subject mask by + # computing the intersection of all run masks. + mask_intersection = Node(MultiImageMaths(), name = 'mask_intersection') + mask_intersection.inputs.op_string = '-mul %s ' * (len(self.run_list) - 1) + subject_level.connect(split_masks, 'out1', mask_intersection, 'in_file') + subject_level.connect(split_masks, 'out2', mask_intersection, 'operand_files') + + # L2Model Node - Generate subject specific second level model + generate_model = Node(L2Model(), name = 'generate_model') + generate_model.inputs.num_copes = len(self.run_list) + + # FLAMEO Node - Estimate model + estimate_model = Node(FLAMEO(), name = 'estimate_model') + estimate_model.inputs.run_mode = 'flame1' + subject_level.connect(mask_intersection, 'out_file', estimate_model, 'mask_file') + subject_level.connect(merge_copes, 'merged_file', estimate_model, 'cope_file') + subject_level.connect(merge_varcopes, 'merged_file', estimate_model, 'var_cope_file') + subject_level.connect(generate_model, 'design_mat', estimate_model, 'design_file') + subject_level.connect(generate_model, 'design_con', estimate_model, 't_con_file') + subject_level.connect(generate_model, 'design_grp', estimate_model, 'cov_split_file') + + # DataSink Node - store the wanted results in the wanted directory + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + subject_level.connect( + mask_intersection, 'out_file', data_sink, 'subject_level_analysis.@mask') + subject_level.connect(estimate_model, 'zstats', data_sink, 'subject_level_analysis.@stats') + subject_level.connect( + estimate_model, 'tstats', data_sink, 'subject_level_analysis.@tstats') + subject_level.connect(estimate_model, 'copes', data_sink, 'subject_level_analysis.@copes') + subject_level.connect( + estimate_model, 'var_copes', data_sink, 'subject_level_analysis.@varcopes') + + return subject_level + + def get_subject_level_outputs(self): + """ Return the names of the files the subject level analysis is supposed to generate. """ + + parameters = { + 'contrast_id' : self.contrast_list, + 'subject_id' : self.subject_list, + 'file' : ['cope1.nii.gz', 'tstat1.nii.gz', 'varcope1.nii.gz', 'zstat1.nii.gz'] + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'subject_level_analysis', '_contrast_id_{contrast_id}_subject_id_{subject_id}','{file}' + ) + return_list = [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + parameters = { + 'contrast_id' : self.contrast_list, + 'subject_id' : self.subject_list, + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'subject_level_analysis', '_contrast_id_{contrast_id}_subject_id_{subject_id}', + 'sub-{subject_id}_task-MGT_run-01_bold_space-MNI152NLin2009cAsym_preproc_brain_mask_maths.nii.gz' + ) + return_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + return return_list + + def get_one_sample_t_test_regressors(subject_list: list) -> dict: + """ + Create dictionary of regressors for one sample t-test group analysis. + + Parameters: + - subject_list: ids of subject in the group for which to do the analysis + + Returns: + - dict containing named lists of regressors. + """ + + return dict(group_mean = [1 for _ in subject_list]) + + def get_two_sample_t_test_regressors( + equal_range_ids: list, + equal_indifference_ids: list, + subject_list: list, + ) -> dict: + """ + Create dictionary of regressors for two sample t-test group analysis. + + Parameters: + - equal_range_ids: ids of subjects in equal range group + - equal_indifference_ids: ids of subjects in equal indifference group + - subject_list: ids of subject for which to do the analysis + + Returns: + - regressors, dict: containing named lists of regressors. + - groups, list: group identifiers to distinguish groups in FSL analysis. + """ + + # Create 2 lists containing n_sub values which are + # * 1 if the participant is on the group + # * 0 otherwise + equal_range_regressors = [1 if i in equal_range_ids else 0 for i in subject_list] + equal_indifference_regressors = [ + 1 if i in equal_indifference_ids else 0 for i in subject_list + ] + + # Create regressors output : a dict with the two list + regressors = dict( + equalRange = equal_range_regressors, + equalIndifference = equal_indifference_regressors + ) + + # Create groups outputs : a list with 1 for equalRange subjects and 2 for equalIndifference + groups = [1 if i == 1 else 2 for i in equal_range_regressors] + + return regressors, groups + + def get_group_level_analysis(self): + """ + Return all workflows for the group level analysis. + + Returns; + - a list of nipype.WorkFlow + """ + + methods = ['equalRange', 'equalIndifference', 'groupComp'] + return [self.get_group_level_analysis_sub_workflow(method) for method in methods] + + def get_group_level_analysis_sub_workflow(self, method): + """ + Return a workflow for the group level analysis. + + Parameters: + - method: one of 'equalRange', 'equalIndifference' or 'groupComp' + + Returns: + - group_level: nipype.WorkFlow + """ + # Compute the number of participants used to do the analysis + nb_subjects = len(self.subject_list) + + # Declare the workflow + group_level = Workflow( + base_dir = self.directories.working_dir, + name = f'group_level_analysis_{method}_nsub_{nb_subjects}') + + # Infosource Node - iterate over the contrasts generated by the subject level analysis + information_source = Node(IdentityInterface( + fields = ['contrast_id']), + name = 'information_source') + information_source.iterables = [('contrast_id', self.contrast_list)] + + # SelectFiles Node - select necessary files + templates = { + 'cope' : join(self.directories.output_dir, + 'subject_level_analysis', '_contrast_id_{contrast_id}_subject_id_*', + 'cope1.nii.gz'), + 'varcope' : join(self.directories.output_dir, + 'subject_level_analysis', '_contrast_id_{contrast_id}_subject_id_*', + 'varcope1.nii.gz'), + 'masks': join(self.directories.output_dir, + 'subject_level_analysis', '_contrast_id_1_subject_id_*', + 'sub-*_task-MGT_run-*_bold_space-MNI152NLin2009cAsym_preproc_brain_mask_maths.nii.gz') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory = self.directories.results_dir + group_level.connect(information_source, 'contrast_id', select_files, 'contrast_id') + + # Function Node elements_in_string + # Get contrast of parameter estimates (cope) for these subjects + # Note : using a MapNode with elements_in_string requires using clean_list to remove + # None values from the out_list + get_copes = MapNode(Function( + function = elements_in_string, + input_names = ['input_str', 'elements'], + output_names = ['out_list'] + ), + name = 'get_copes', iterfield = 'input_str' + ) + group_level.connect(select_files, 'cope', get_copes, 'input_str') + + # Function Node elements_in_string + # Get variance of the estimated copes (varcope) for these subjects + # Note : using a MapNode with elements_in_string requires using clean_list to remove + # None values from the out_list + get_varcopes = MapNode(Function( + function = elements_in_string, + input_names = ['input_str', 'elements'], + output_names = ['out_list'] + ), + name = 'get_varcopes', iterfield = 'input_str' + ) + group_level.connect(select_files, 'varcope', get_varcopes, 'input_str') + + # Merge Node - Merge cope files + merge_copes = Node(Merge(), name = 'merge_copes') + merge_copes.inputs.dimension = 't' + group_level.connect(get_copes, ('out_list', clean_list), merge_copes, 'in_files') + + # Merge Node - Merge cope files + merge_varcopes = Node(Merge(), name = 'merge_varcopes') + merge_varcopes.inputs.dimension = 't' + group_level.connect(get_varcopes, ('out_list', clean_list), merge_varcopes, 'in_files') + + # Split Node - Split mask list to serve them as inputs of the MultiImageMaths node. + split_masks = Node(Split(), name = 'split_masks') + split_masks.inputs.splits = [1, len(self.subject_list) - 1] + split_masks.inputs.squeeze = True # Unfold one-element splits removing the list + group_level.connect(select_files, 'masks', split_masks, 'inlist') + + # MultiImageMaths Node - Create a subject mask by + # computing the intersection of all run masks. + mask_intersection = Node(MultiImageMaths(), name = 'mask_intersection') + mask_intersection.inputs.op_string = '-mul %s ' * (len(self.subject_list) - 1) + group_level.connect(split_masks, 'out1', mask_intersection, 'in_file') + group_level.connect(split_masks, 'out2', mask_intersection, 'operand_files') + + # MultipleRegressDesign Node - Specify model + specify_model = Node(MultipleRegressDesign(), name = 'specify_model') + + # FLAMEO Node - Estimate model + estimate_model = Node(FLAMEO(), name = 'estimate_model') + estimate_model.inputs.run_mode = 'flame1' + group_level.connect(mask_intersection, 'out_file', estimate_model, 'mask_file') + group_level.connect(merge_copes, estimate_model, [('merged_file', 'cope_file')]), + group_level.connect(merge_varcopes, estimate_model, [('merged_file', 'var_cope_file')]), + group_level.connect(specify_model, 'design_mat', estimate_model, 'design_file') + group_level.connect(specify_model, 'design_con', estimate_model, 't_con_file') + group_level.connect(specify_model, 'design_grp', estimate_model, 'cov_split_file') + + # Randomise Node - Perform clustering on statistical output + randomise = MapNode(Randomise(), + name = 'randomise', + iterfield = ['in_file', 'dlh', 'volume', 'cope_file'], + synchronize = True) + randomise.inputs.tfce = True + randomise.inputs.num_perm = 10000 + randomise.inputs. = True + randomise.inputs.pthreshold = 0.05 + + +tcon +mask +fcon +design_mat + + (estimate_model, randomise, + ('zstats', 'in_file') + ('copes', 'cope_file')) + + + + + # Datasink Node - save important files + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + (estimate_model, data_sink, [ + ('zstats', f'group_level_analysis_{method}_nsub_{nb_subjects}.@zstats'), + ('tstats', f'group_level_analysis_{method}_nsub_{nb_subjects}.@tstats') + ]), + (cluster, data_sink, [ + ('threshold_file', f'group_level_analysis_{method}_nsub_{nb_subjects}.@thresh'), + ('pval_file', f'group_level_analysis_{method}_nsub_{nb_subjects}.@pval')]) + ]) + + if method in ('equalIndifference', 'equalRange'): + # Setup a one sample t-test + specify_model.inputs.contrasts = [ + ['group_mean', 'T', ['group_mean'], [1]], + ['group_mean_neg', 'T', ['group_mean'], [-1]] + ] + + # Function Node get_group_subjects - Get subjects in the group and in the subject_list + get_group_subjects = Node(Function( + function = list_intersection, + input_names = ['list_1', 'list_2'], + output_names = ['out_list'] + ), + name = 'get_group_subjects' + ) + get_group_subjects.inputs.list_1 = get_group(method) + get_group_subjects.inputs.list_2 = self.subject_list + + # Function Node get_one_sample_t_test_regressors + # Get regressors in the equalRange and equalIndifference method case + regressors_one_sample = Node( + Function( + function = self.get_one_sample_t_test_regressors, + input_names = ['subject_list'], + output_names = ['regressors'] + ), + name = 'regressors_one_sample', + ) + + # Add missing connections + group_level_analysis.connect([ + (get_group_subjects, get_copes, [('out_list', 'elements')]), + (get_group_subjects, get_varcopes, [('out_list', 'elements')]), + (get_group_subjects, regressors_one_sample, [('out_list', 'subject_list')]), + (regressors_one_sample, specify_model, [('regressors', 'regressors')]) + ]) + + elif method == 'groupComp': + + # Select copes and varcopes corresponding to the selected subjects + # Indeed the SelectFiles node asks for all (*) subjects available + get_copes.inputs.elements = self.subject_list + get_varcopes.inputs.elements = self.subject_list + + # Setup a two sample t-test + specify_model.inputs.contrasts = [ + ['equalRange_sup', 'T', ['equalRange', 'equalIndifference'], [1, -1]] + ] + + # Function Node get_equal_range_subjects + # Get subjects in the equalRange group and in the subject_list + get_equal_range_subjects = Node(Function( + function = list_intersection, + input_names = ['list_1', 'list_2'], + output_names = ['out_list'] + ), + name = 'get_equal_range_subjects' + ) + get_equal_range_subjects.inputs.list_1 = get_group('equalRange') + get_equal_range_subjects.inputs.list_2 = self.subject_list + + # Function Node get_equal_indifference_subjects + # Get subjects in the equalIndifference group and in the subject_list + get_equal_indifference_subjects = Node(Function( + function = list_intersection, + input_names = ['list_1', 'list_2'], + output_names = ['out_list'] + ), + name = 'get_equal_indifference_subjects' + ) + get_equal_indifference_subjects.inputs.list_1 = get_group('equalIndifference') + get_equal_indifference_subjects.inputs.list_2 = self.subject_list + + # Function Node get_two_sample_t_test_regressors + # Get regressors in the groupComp method case + regressors_two_sample = Node( + Function( + function = self.get_two_sample_t_test_regressors, + input_names = [ + 'equal_range_ids', + 'equal_indifference_ids', + 'subject_list', + ], + output_names = ['regressors', 'groups'] + ), + name = 'regressors_two_sample', + ) + regressors_two_sample.inputs.subject_list = self.subject_list + + # Add missing connections + group_level_analysis.connect([ + (get_equal_range_subjects, regressors_two_sample, [ + ('out_list', 'equal_range_ids') + ]), + (get_equal_indifference_subjects, regressors_two_sample, [ + ('out_list', 'equal_indifference_ids') + ]), + (regressors_two_sample, specify_model, [ + ('regressors', 'regressors'), + ('groups', 'groups')]) + ]) + + return group_level_analysis + + def get_group_level_outputs(self): + """ Return all names for the files the group level analysis is supposed to generate. """ + + # Handle equalRange and equalIndifference + parameters = { + 'contrast_id': self.contrast_list, + 'method': ['equalRange', 'equalIndifference'], + 'file': [ + '_cluster0/zstat1_pval.nii.gz', + '_cluster0/zstat1_threshold.nii.gz', + '_cluster1/zstat2_pval.nii.gz', + '_cluster1/zstat2_threshold.nii.gz', + 'tstat1.nii.gz', + 'tstat2.nii.gz', + 'zstat1.nii.gz', + 'zstat2.nii.gz' + ] + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'group_level_analysis_{method}_nsub_'+f'{len(self.subject_list)}', + '_contrast_id_{contrast_id}', + '{file}' + ) + return_list = [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + # Handle groupComp + parameters = { + 'contrast_id': self.contrast_list, + 'file': [ + '_cluster0/zstat1_pval.nii.gz', + '_cluster0/zstat1_threshold.nii.gz', + 'tstat1.nii.gz', + 'zstat1.nii.gz' + ] + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + f'group_level_analysis_groupComp_nsub_{len(self.subject_list)}', + '_contrast_id_{contrast_id}', + '{file}' + ) + return_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + return return_list + + def get_hypotheses_outputs(self): + """ Return all hypotheses output file names. """ + + nb_sub = len(self.subject_list) + files = [ + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz'), + join(f'group_level_analysis_groupComp_nsub_{nb_sub}', + '_contrast_id_2', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_groupComp_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz') + ] + return [join(self.directories.output_dir, f) for f in files] diff --git a/tests/pipelines/test_team_B5I6.py b/tests/pipelines/test_team_B5I6.py new file mode 100644 index 00000000..f0ed393f --- /dev/null +++ b/tests/pipelines/test_team_B5I6.py @@ -0,0 +1,187 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Tests of the 'narps_open.pipelines.team_B5I6' module. + +Launch this test with PyTest + +Usage: +====== + pytest -q test_team_B5I6.py + pytest -q test_team_B5I6.py -k +""" +from os.path import join, exists, abspath +from filecmp import cmp + +from pytest import helpers, mark +from nipype import Workflow, Node, Function +from nipype.interfaces.base import Bunch + +from narps_open.utils.configuration import Configuration +from narps_open.pipelines.team_B5I6 import PipelineTeamB5I6 + +class TestPipelinesTeamB5I6: + """ A class that contains all the unit tests for the PipelineTeamB5I6 class.""" + + @staticmethod + @mark.unit_test + def test_create(): + """ Test the creation of a PipelineTeamB5I6 object """ + + pipeline = PipelineTeamB5I6() + + # 1 - check the parameters + assert pipeline.fwhm == 5.0 + assert pipeline.team_id == 'B5I6' + + # 2 - check workflows + assert pipeline.get_preprocessing() is None + assert isinstance(pipeline.get_run_level_analysis(), Workflow) + assert isinstance(pipeline.get_subject_level_analysis(), Workflow) + group_level = pipeline.get_group_level_analysis() + assert len(group_level) == 3 + for sub_workflow in group_level: + assert isinstance(sub_workflow, Workflow) + + @staticmethod + @mark.unit_test + def test_outputs(): + """ Test the expected outputs of a PipelineTeamB5I6 object """ + + pipeline = PipelineTeamB5I6() + + # 1 - 1 subject outputs + pipeline.subject_list = ['001'] + helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*1, 4*4*1 + 4*1, 8*4*2 + 4*4, 18]) + + # 2 - 4 subjects outputs + pipeline.subject_list = ['001', '002', '003', '004'] + helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*4, 4*4*4 + 4*4, 8*4*2 + 4*4, 18]) + + @staticmethod + @mark.unit_test + def test_subject_information(): + """ Test the get_subject_information method """ + + # Get test files + test_file = join(Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') + test_file_2 = join( + Configuration()['directories']['test_data'], 'pipelines', 'events_resp.tsv') + + # Prepare several scenarii + info_missed = PipelineTeamB5I6.get_subject_information(test_file) + info_no_missed = PipelineTeamB5I6.get_subject_information(test_file_2) + + # Compare bunches to expected + bunch = info_missed[0] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['trial', 'missed'] + helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 27.535, 36.435], [19.535]]) + helpers.compare_float_2d_arrays(bunch.durations, [[4.0, 4.0, 4.0, 4.0], [4.0]]) + assert bunch.amplitudes is None + assert bunch.tmod is None + assert bunch.regressor_names is None + assert bunch.regressors is None + bunch = bunch.pmod[0] + assert isinstance(bunch, Bunch) + assert bunch.name == ['gain', 'loss'] + assert bunch.poly == [1, 1] + helpers.compare_float_2d_arrays(bunch.param, [ + [-4.5, 15.5, -8.5, -2.5], + [-7.0, 1.0, 2.0, 4.0] + ]) + + bunch = info_no_missed[0] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['trial'] + helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 27.535, 36.435]]) + helpers.compare_float_2d_arrays(bunch.durations, [[4.0, 4.0, 4.0, 4.0]]) + assert bunch.amplitudes is None + assert bunch.tmod is None + assert bunch.regressor_names is None + assert bunch.regressors is None + bunch = bunch.pmod[0] + assert isinstance(bunch, Bunch) + assert bunch.name == ['gain', 'loss'] + assert bunch.poly == [1, 1] + helpers.compare_float_2d_arrays(bunch.param, [ + [-4.5, 15.5, -8.5, -2.5], + [-7.0, 1.0, 2.0, 4.0] + ]) + + @staticmethod + @mark.unit_test + def test_confounds_file_no_outliers(temporary_data_dir): + """ Test the get_confounds_file method in the case with no outliers """ + + # Get input and reference output file + confounds_file = join( + Configuration()['directories']['test_data'], 'pipelines', 'confounds.tsv') + reference_file = join( + Configuration()['directories']['test_data'], + 'pipelines', 'team_B5I6', 'out_confounds_no_outliers.tsv') + + # Create new confounds file + confounds_node = Node(Function( + input_names = ['filepath', 'subject_id', 'run_id'], + output_names = ['confounds_file'], + function = PipelineTeamB5I6.get_confounds_file), + name = 'confounds_node') + confounds_node.base_dir = temporary_data_dir + confounds_node.inputs.filepath = confounds_file + confounds_node.inputs.subject_id = 'sid' + confounds_node.inputs.run_id = 'rid' + confounds_node.run() + + # Check confounds file was created + created_confounds_file = abspath(join( + temporary_data_dir, confounds_node.name, 'confounds_file_sub-sid_run-rid.tsv')) + assert exists(created_confounds_file) + + from shutil import copy + copy(created_confounds_file, '/work/toot.tsv') + + # Check contents + assert cmp(reference_file, created_confounds_file) + + @staticmethod + @mark.unit_test + def test_confounds_file_outliers(temporary_data_dir): + """ Test the get_confounds_file method in the case with outliers """ + + # Get input and reference output file + confounds_file = join( + Configuration()['directories']['test_data'], + 'pipelines', 'team_B5I6', 'confounds_with_outliers.tsv') + reference_file = join( + Configuration()['directories']['test_data'], + 'pipelines', 'team_B5I6', 'out_confounds_outliers.tsv') + + # Create new confounds file + confounds_node = Node(Function( + input_names = ['filepath', 'subject_id', 'run_id'], + output_names = ['confounds_file'], + function = PipelineTeamB5I6.get_confounds_file), + name = 'confounds_node') + confounds_node.base_dir = temporary_data_dir + confounds_node.inputs.filepath = confounds_file + confounds_node.inputs.subject_id = 'sid' + confounds_node.inputs.run_id = 'rid' + confounds_node.run() + + # Check confounds file was created + created_confounds_file = abspath(join( + temporary_data_dir, confounds_node.name, 'confounds_file_sub-sid_run-rid.tsv')) + assert exists(created_confounds_file) + + from shutil import copy + copy(created_confounds_file, '/work/toot_outl.tsv') + + # Check contents + assert cmp(reference_file, created_confounds_file) + + @staticmethod + @mark.pipeline_test + def test_execution(): + """ Test the execution of a PipelineTeamB5I6 and compare results """ + helpers.test_pipeline_evaluation('B5I6') diff --git a/tests/test_data/pipelines/team_B5I6/confounds_with_outliers.tsv b/tests/test_data/pipelines/team_B5I6/confounds_with_outliers.tsv new file mode 100644 index 00000000..5ff676fe --- /dev/null +++ b/tests/test_data/pipelines/team_B5I6/confounds_with_outliers.tsv @@ -0,0 +1,20 @@ +CSF WhiteMatter GlobalSignal stdDVARS non-stdDVARS vx-wisestdDVARS FramewiseDisplacement tCompCor00 tCompCor01 tCompCor02 tCompCor03 tCompCor04 tCompCor05 aCompCor00 aCompCor01 aCompCor02 aCompCor03 aCompCor04 aCompCor05 Cosine00 Cosine01 Cosine02 Cosine03 Cosine04 Cosine05 NonSteadyStateOutlier00 X Y Z RotX RotY RotZ +6551.281999999999 6476.4653 9874.576 n/a n/a n/a n/a 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 -0.0 0.0 +6484.7285 6473.4890000000005 9830.212 1.09046686 52.78273392 1.05943739 0.13527900930999998 0.0263099209 -0.0673065879 0.0934882554 -0.0079328884 0.0338007737 -0.011491083999999999 -0.042411347099999996 0.027736422900000002 0.0453303087 -0.07022609490000001 0.0963618709 -0.0200867957 0.0665186088 0.0665174038 0.0665153954 0.0665125838 0.0665089688 0.06650455059999999 0.0 -0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 +6441.5337 6485.7256 9821.212 1.07520139 152.04382706 1.03821933 0.12437666391 -0.0404820317 0.034150583 0.13661184210000002 0.0745358691 -0.0054829985999999995 -0.0217322686 0.046214115199999996 0.005774624 -0.043909359800000006 -0.075619539 0.17546891539999998 -0.0345256763 0.0665153954 0.06650455059999999 0.06648647719999999 0.0664611772 0.0664286533 0.0663889091 0.0 -2.56954e-05 -0.00923735 0.0549667 0.000997278 -0.00019745 -0.000398988 +6401.115 6481.4287 9808.847 1.05443561 51.03868103 1.024683 0.14159200100000002 0.0416039596 0.043600489299999996 0.0194078842 -0.0584053666 -0.0393477511 -0.0862903833 0.0048847126000000005 0.0524381164 -0.0095181227 -0.026580354100000002 0.045246431 0.0685650559 0.0665089688 0.0664788467 0.0664286533 0.0663584038 0.0662681193 0.0661578271 0.0 -3.28464e-05 0.0579411 0.0692827 0.000114534 -3.75912e-24 -0.000277374 +6428.2754 6486.5215 9800.982 1.03301764 50.0019722 1.00801814 0.0500087547 0.039141912 -0.048557309900000005 0.0023491611 0.0570407992 0.056129457800000004 -0.0848882627 -0.0639463374 0.006067867199999999 -0.005321996 0.1230282776 0.0411739245 0.11016770869999999 0.0664993292 0.0664402971 0.0663419489 0.0662043429 0.0660275606 0.0658117066 0.0 -2.82417e-05 0.0666809 0.0471475 -0.000130534 -0.0 -0.000139859 +6379.1426 6479.1035 9790.05 1.067186 51.65584946 1.0255338 0.1444939117 -0.0570486116 0.0139479967 0.0350964936 0.0615363225 -0.0546076932 -0.0216442783 0.040404627 -0.0172438568 -0.0327785129 0.0366254628 0.0290480775 -0.012212245200000001 0.06648647719999999 0.0663889091 0.0662264016 0.0659991137 0.0657072678 0.0653511492 0.0 0.00212317 0.00537855 0.0604605 0.000938976 -5.98122e-23 -0.000424892 +6389.1616 6485.956999999999 9798.146999999999 1.03589809 50.14139557 1.01152003 0.06430312549999999 0.0488044693 0.04843859 0.019325393 -0.051024495499999996 -0.0774132526 -0.0439173899 -0.0135618257 0.0278786709 -0.0400976257 -0.0845411567 -0.0046610891 0.0438297914 0.0664704133 0.0663246928 0.0660820618 0.0657428748 0.0653076276 0.0647769559 0.0 -8.45055e-05 0.0420387 0.0650003 0.000522193 -4.5944699999999996e-23 -0.000423765 +6416.148 6487.11 9805.799 1.03183794 49.94487 1.00848448 0.0877196945 0.046519325099999995 -0.045735816799999995 0.0531509041 0.033536665099999995 0.0426497861 -0.0671152129 -0.0687792575 0.006283962 -0.0144780862 -0.0259729938 -0.0291949146 0.014316824799999999 0.0664511384 0.0662476605 0.0659089921 0.06543582419999999 0.0648291226 0.0640901252 0.0 -0.0219226 0.0512385 0.0550966 0.000272193 -0.00046169 -0.000199893 +6407.106 6487.23 9806.202 1.06404078 51.50361252 1.02338874 0.12704458999999999 -0.0628035146 -0.0135729873 0.034224742999999995 0.0181556141 -0.0597216613 -0.0524812854 0.0377449599 -0.0579724226 -0.0045264773 -0.0370587968 -0.0343268544 -0.079852853 0.0664286533 0.0661578271 0.0657072678 0.0650781993 0.0642723306 0.0632918512 0.0 -0.0218915 -0.00836779 0.0631127 0.0012380300000000004 -0.000439598 0.0 +6381.904 6491.920999999999 9801.41 1.04746306 50.70118332 1.0135237 0.07707824 0.0163919608 0.0696230578 -0.0045920252000000005 -0.0363335399 0.043117322 -0.09109143060000001 0.038139249199999996 0.0604562842 -0.0755071391 0.006743919 -0.0060393816 -0.0042580964 0.0664029592 0.0660552101 0.0654769766 0.0646702763 0.0636379241 0.0623835221 0.0 -0.0218846 0.0451283 0.0698908 0.0009294180000000002 -0.000466929 0.0 +6443.9834 6473.995 9794.653 1.03327584 50.01446915 1.00756884 0.033159035 0.037253001699999996 -0.0254029516 -0.0782083535 0.0647054328 0.0508582046 -0.0976487128 -0.034160880899999996 -0.0040617346 0.0487996336 0.0848041487 -0.0018666815 0.055534395300000004 0.0663740572 0.0659398291 0.0652182186 0.06421237070000001 0.0629266692 0.0613667174 0.0 -0.0218922 0.0389939 0.0499376 0.000936915 -0.00053235 -6.83587e-05 +6399.231 6488.91 9777.269 1.0600059 51.30830383 1.02645695 0.124505895 -0.0608352655 -0.038347576099999996 -0.0783155209 -0.047988137300000004 -0.012615503600000001 -0.0042425137 0.0181519689 -0.048044653 -0.0082978382 0.013016688 -0.1543962652 -0.0265575586 0.0663419489 0.0658117066 0.0649311064 0.0637048363 0.0621394246 0.060243205 0.0 -0.0218563 -0.00877046 0.0643078 0.00199082 -0.000656795 0.0 +6355.352 6477.0786 9755.872 1.06057775 51.33598709 1.02456272 0.08770730999999998 -0.0080858687 0.0595016954 -0.0490129325 -0.0124959526 0.0682339906 0.0154944714 0.0590289672 0.0287994709 0.0346818337 -0.010791274299999999 0.009191705 0.055267397 0.0663066357 0.0656708672 0.0646157647 0.0631480654 0.0612771412 0.0590149387 0.0 -0.0218688 0.0371576 0.0722054 0.00146155 -0.000508682 0.0 +6356.6245 6467.55 9748.421 1.02207148 49.47213745 1.00191152 0.0205794 0.0385542551 -0.0346451025 -0.0747522904 0.060552027599999995 0.0626499417 0.0415792814 -0.0460955165 0.0104297412 0.10550639560000001 0.0556434251 -0.044798088 0.0105382343 0.0662681193 0.0655173382 0.0642723306 0.0625424883 0.0603408602 0.0576840542 0.0 -0.0218522 0.0319276 0.0612045 0.00146155 -0.00059532 0.0 +6326.0557 6465.569 9738.453000000001 1.05498719 51.06538391 1.01648581 0.07538484 -0.075319679 0.0006298104 -0.032952762999999996 0.044291419299999996 -0.049490735099999995 0.080840249 0.0595743027 -0.028879615 0.0431988675 -0.0414852927 -0.079980074 -0.09799051019999999 0.0662264016 0.0653511492 0.0639009535 0.0618885731 0.0593317122 0.0562528657 0.0 -0.0318116 0.00937896 0.0779885 0.00192327 -0.000655456 0.0 +6374.3335 6472.7896 9753.27 1.04270554 50.47090149 1.01588619 0.10981649 0.059109603 0.0398218045 0.0082860844 -0.054184145700000005 -0.0681863947 0.0432304277 -0.007844685 0.0578450393 0.0055081056 0.0113240728 0.0268679277 -0.031792571 0.0661814847 0.0651723324 0.0635017949 0.061186825099999995 0.0582509159 0.054723861799999995 0.0 -0.0276812 0.0703663 0.084403 0.00133038 -0.000635065 -0.000152404 +6378.758000000001 6469.625 9767.294 1.02955639 49.83443451 1.00836587 0.0392554 -0.0270312077 -0.0381127912 0.0381207645 -0.0829419772 0.006719891800000001 -0.0079905833 -0.047767327000000005 -0.0703481828 -0.0014303484 -0.0998887326 -0.0148952053 -0.08504836060000001 0.06613337059999999 0.06498092230000001 0.0630750281 0.0604377867 0.057099776500000005 0.053099701299999996 0.0 -0.0277335 0.0349253 0.08561299999999998 0.00133038 -0.000656881 -0.000123178 +6371.084 6468.0317 9774.892 1.03255022 49.97934723 1.00839365 0.02700655000000001 -0.0020486699 0.0808370966 0.035830127 -0.0338654606 0.039639693399999995 -0.0005207792 0.0357297151 0.0434899465 -0.0110042022 -0.1272102049 -0.0595300301 -0.0299245315 0.0660820618 0.0647769559 0.0626208389 0.0596420368 0.0558796839 0.051383208300000004 0.0 -0.0295774 0.0316995 0.0924462 0.00145271 -0.000600316 0.0 +6400.2603 6480.2183 9781.016 1.02676523 49.69933319 1.00593698 10.1047363 0.032433138 -0.029600444500000003 -0.0303759073 0.043935624000000006 0.0541548529 -0.026322396600000002 -0.0382086532 0.013210211299999998 -0.0063683262 0.048276844000000006 0.0736115089 0.052873248 0.0660275606 0.0645604726 0.0621394246 0.058800190599999996 0.0545921117 0.0495773675 0.0 -0.0281078 0.0837567 0.0793704 0.000708306 -0.000618586 0.0 diff --git a/tests/test_data/pipelines/team_B5I6/out_confounds_no_outliers.tsv b/tests/test_data/pipelines/team_B5I6/out_confounds_no_outliers.tsv new file mode 100644 index 00000000..1453081d --- /dev/null +++ b/tests/test_data/pipelines/team_B5I6/out_confounds_no_outliers.tsv @@ -0,0 +1,3 @@ +0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +-0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 -0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 9.937996410250001e-05 0.00098247141136 9.0559466761e-12 1.7605839969e-06 1.4760426124900002e-07 2.82878761e-08 9.937996410250001e-05 0.00098247141136 9.0559466761e-12 1.7605839969e-06 1.4760426124900002e-07 2.82878761e-08 +-2.56954e-05 -0.00923735 0.0549667 0.000997278 -0.00019745 -0.000398988 0.009943254600000001 0.022107050000000003 0.05496970931 -0.00032959199999999986 0.000186743 -0.00023079800000000002 6.6025358116e-10 8.53286350225e-05 0.00302133810889 9.94563409284e-07 3.89865025e-08 1.59191424144e-07 9.886831204042119e-05 0.0004887216597025001 0.003021668941625901 1.0863088646399991e-07 3.4872948049e-08 5.326771680400001e-08 diff --git a/tests/test_data/pipelines/team_B5I6/out_confounds_outliers.tsv b/tests/test_data/pipelines/team_B5I6/out_confounds_outliers.tsv new file mode 100644 index 00000000..face8e1b --- /dev/null +++ b/tests/test_data/pipelines/team_B5I6/out_confounds_outliers.tsv @@ -0,0 +1,19 @@ +0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +-0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 -0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 9.937996410250001e-05 0.00098247141136 9.0559466761e-12 1.7605839969e-06 1.4760426124900002e-07 2.82878761e-08 9.937996410250001e-05 0.00098247141136 9.0559466761e-12 1.7605839969e-06 1.4760426124900002e-07 2.82878761e-08 0.0 +-2.56954e-05 -0.00923735 0.0549667 0.000997278 -0.00019745 -0.000398988 0.009943254600000001 0.022107050000000003 0.05496970931 -0.00032959199999999986 0.000186743 -0.00023079800000000002 6.6025358116e-10 8.53286350225e-05 0.00302133810889 9.94563409284e-07 3.89865025e-08 1.59191424144e-07 9.886831204042119e-05 0.0004887216597025001 0.003021668941625901 1.0863088646399991e-07 3.4872948049e-08 5.326771680400001e-08 1.0 +-3.28464e-05 0.0579411 0.0692827 0.000114534 -3.75912e-24 -0.000277374 -7.150999999999996e-06 0.06717845 0.014316000000000002 -0.0008827440000000001 0.00019745 0.00012161400000000003 1.0788859929599998e-09 0.00335717106921 0.0048000925192900005 1.3118037155999999e-08 1.41309831744e-47 7.6936335876e-08 5.113680099999995e-11 0.0045129441444025 0.00020494785600000008 7.792369695360002e-07 3.89865025e-08 1.4789964996000009e-08 0.0 +-2.82417e-05 0.0666809 0.0471475 -0.000130534 -0.0 -0.000139859 4.604699999999997e-06 0.008739799999999999 -0.0221352 -0.000245068 3.75912e-24 0.000137515 7.975936188900001e-10 0.00444634242481 0.00222288675625 1.7039125155999998e-08 0.0 1.9560539880999998e-08 2.120326208999997e-11 7.638410403999998e-05 0.00048996707904 6.005832462399999e-08 1.41309831744e-47 1.8910375225e-08 0.0 +0.00212317 0.00537855 0.0604605 0.000938976 -5.98122e-23 -0.000424892 0.0021514117 -0.06130235 0.013312999999999998 0.00106951 -5.98122e-23 -0.000285033 4.5078508489000005e-06 2.8928800102500002e-05 0.00365547206025 8.81675928576e-07 3.5774992688400005e-45 1.80533211664e-07 4.62857230289689e-06 0.0037579781155225 0.00017723596899999996 1.1438516401000002e-06 3.5774992688400005e-45 8.124381108900001e-08 0.0 +-8.45055e-05 0.0420387 0.0650003 0.000522193 -4.59447e-23 -0.000423765 -0.0022076755000000003 0.036660149999999996 0.004539799999999997 -0.000416783 1.38675e-23 1.1270000000000182e-06 7.14117953025e-09 0.0017672522976899998 0.00422503900009 2.72685529249e-07 2.1109154580900002e-45 1.7957677522499998e-07 4.8738311133002515e-06 0.0013439665980224996 2.060978403999997e-05 1.73708069089e-07 1.9230755625e-46 1.270129000000041e-12 0.0 +-0.0219226 0.0512385 0.0550966 0.000272193 -0.00046169 -0.000199893 -0.021838094500000002 0.009199800000000001 -0.009903699999999994 -0.00025 -0.00046169 0.000223872 0.00048060039076 0.00262538388225 0.0030356353315600004 7.4089029249e-08 2.131576561e-07 3.9957211449000005e-08 0.00047690237139093035 8.463632004000002e-05 9.808327368999989e-05 6.25e-08 2.131576561e-07 5.011867238399999e-08 0.0 +-0.0218915 -0.00836779 0.0631127 0.00123803 -0.000439598 0.0 3.1099999999999184e-05 -0.05960629 0.008016099999999991 0.000965837 2.2092000000000032e-05 0.000199893 0.00047923777225000004 7.00199094841e-05 0.003983212901289999 1.5327182809e-06 1.9324640160399999e-07 0.0 9.672099999999491e-10 0.0035529098075640997 6.425785920999987e-05 9.32841110569e-07 4.880564640000014e-10 3.9957211449000005e-08 0.0 +-0.0218846 0.0451283 0.0698908 0.000929418 -0.000466929 0.0 6.900000000000656e-06 0.05349609 0.006778100000000009 -0.000308612 -2.7331000000000013e-05 0.0 0.00047893571716 0.00203656346089 0.0048847239246400005 8.63817818724e-07 2.1802269104099998e-07 0.0 4.7610000000009045e-11 0.0028618316452881003 4.5942639610000124e-05 9.524136654399999e-08 7.469835610000007e-10 0.0 0.0 +-0.0218922 0.0389939 0.0499376 0.000936915 -0.00053235 -6.83587e-05 -7.599999999999968e-06 -0.006134400000000005 -0.019953200000000004 7.4969999999999724e-06 -6.542100000000002e-05 -6.83587e-05 0.00047926842084000003 0.0015205242372099998 0.00249376389376 8.77809717225e-07 2.833965225e-07 4.672911865690001e-09 5.775999999999951e-11 3.7630863360000066e-05 0.0003981301902400002 5.6205008999999586e-11 4.279907241000003e-09 4.672911865690001e-09 0.0 +-0.0218563 -0.00877046 0.0643078 0.00199082 -0.000656795 0.0 3.59000000000019e-05 -0.04776436 0.0143702 0.001053905 -0.00012444499999999994 6.83587e-05 0.00047769784968999996 7.692096861160001e-05 0.00413549314084 3.9633642723999995e-06 4.3137967202499993e-07 0.0 1.2888100000001365e-09 0.0022814340862096 0.00020650264803999998 1.1107157490249998e-06 1.5486558024999985e-08 4.672911865690001e-09 0.0 +-0.0218688 0.0371576 0.0722054 0.00146155 -0.000508682 0.0 -1.2500000000002093e-05 0.04592806 0.007897600000000005 -0.0005292699999999999 0.00014811299999999996 0.0 0.00047824441344 0.0013806872377599999 0.0052136197891600004 2.1361284025000002e-06 2.58757377124e-07 0.0 1.5625000000005232e-10 0.0021093866953636 6.237208576000007e-05 2.801267328999999e-07 2.1937460768999988e-08 0.0 0.0 +-0.0218522 0.0319276 0.0612045 0.00146155 -0.00059532 0.0 1.660000000000203e-05 -0.0052299999999999985 -0.011000900000000001 0.0 -8.663799999999999e-05 0.0 0.00047751864483999994 0.00101937164176 0.00374599082025 2.1361284025000002e-06 3.544059024e-07 0.0 2.755600000000674e-10 2.7352899999999983e-05 0.00012101980081000002 0.0 7.506143043999999e-09 0.0 0.0 +-0.0318116 0.00937896 0.0779885 0.00192327 -0.000655456 0.0 -0.009959400000000004 -0.022548640000000002 0.016784 0.0004617199999999999 -6.013600000000002e-05 0.0 0.0010119778945600001 8.796489068160001e-05 0.006082206132250001 3.6989674929e-06 4.29622567936e-07 0.0 9.918964836000008e-05 0.0005084411658496001 0.00028170265600000003 2.131853583999999e-07 3.6163384960000026e-09 0.0 0.0 +-0.0276812 0.0703663 0.084403 0.00133038 -0.000635065 -0.000152404 0.004130400000000003 0.06098734000000001 0.0064145000000000035 -0.00059289 2.039100000000003e-05 -0.000152404 0.0007662488334399999 0.004951416175690001 0.007123866409000001 1.7699109443999997e-06 4.0330755422499995e-07 2.3226979216e-08 1.706020416000002e-05 0.0037194556402756008 4.114581025000004e-05 3.5151855209999997e-07 4.1579288100000116e-10 2.3226979216e-08 0.0 +-0.0277335 0.0349253 0.0856129999999999 0.00133038 -0.000656881 -0.000123178 -5.230000000000165e-05 -0.03544100000000001 0.0012099999999998917 0.0 -2.1816000000000023e-05 2.9226000000000008e-05 0.0007691470222500001 0.00121977658009 0.007329585768999982 1.7699109443999997e-06 4.31492648161e-07 1.5172819683999998e-08 2.735290000000173e-09 0.0012560644810000006 1.464099999999738e-06 0.0 4.759378560000009e-10 8.541590760000004e-10 0.0 +-0.0295774 0.0316995 0.0924462 0.00145271 -0.000600316 0.0 -0.001843899999999999 -0.003225800000000001 0.006833200000000109 0.00012233000000000014 5.656499999999996e-05 0.000123178 0.00087482259076 0.00100485830025 0.008546299894440001 2.1103663441000003e-06 3.6037929985600006e-07 0.0 3.3999672099999965e-06 1.0405785640000006e-05 4.669262224000149e-05 1.4964628900000034e-08 3.199599224999995e-09 1.5172819683999998e-08 0.0 +-0.0281078 0.0837567 0.0793704 0.000708306 -0.000618586 0.0 0.0014696000000000015 0.052057200000000005 -0.013075800000000012 -0.000744404 -1.8269999999999918e-05 0.0 0.0007900484208399999 0.007015184794890001 0.006299660396159999 5.01697389636e-07 3.8264863939599996e-07 0.0 2.1597241600000045e-06 0.0027099520718400004 0.00017097654564000032 5.54137315216e-07 3.33792899999997e-10 0.0 1.0 From 54d04b8f7781fb9e610e1f5014d84912b2ce7caf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Thu, 11 Apr 2024 17:22:45 +0200 Subject: [PATCH 04/12] Completed first version of the pipeline --- narps_open/pipelines/team_B5I6.py | 82 ++++++++++++------------------- tests/pipelines/test_team_B5I6.py | 10 +--- 2 files changed, 34 insertions(+), 58 deletions(-) diff --git a/narps_open/pipelines/team_B5I6.py b/narps_open/pipelines/team_B5I6.py index ad201a9c..bdc65f3a 100644 --- a/narps_open/pipelines/team_B5I6.py +++ b/narps_open/pipelines/team_B5I6.py @@ -12,7 +12,7 @@ from nipype.interfaces.fsl import ( IsotropicSmooth, Level1Design, FEATModel, L2Model, Merge, FLAMEO, FILMGLS, MultipleRegressDesign, - Cluster, SmoothEstimate, FSLCommand, Randomise + FSLCommand, Randomise ) from nipype.algorithms.modelgen import SpecifyModel from nipype.interfaces.fsl.maths import MultiImageMaths @@ -595,8 +595,8 @@ def get_group_level_analysis_sub_workflow(self, method): estimate_model = Node(FLAMEO(), name = 'estimate_model') estimate_model.inputs.run_mode = 'flame1' group_level.connect(mask_intersection, 'out_file', estimate_model, 'mask_file') - group_level.connect(merge_copes, estimate_model, [('merged_file', 'cope_file')]), - group_level.connect(merge_varcopes, estimate_model, [('merged_file', 'var_cope_file')]), + group_level.connect(merge_copes, 'merged_file', estimate_model, 'cope_file') + group_level.connect(merge_varcopes, 'merged_file', estimate_model, 'var_cope_file') group_level.connect(specify_model, 'design_mat', estimate_model, 'design_file') group_level.connect(specify_model, 'design_con', estimate_model, 't_con_file') group_level.connect(specify_model, 'design_grp', estimate_model, 'cov_split_file') @@ -608,33 +608,23 @@ def get_group_level_analysis_sub_workflow(self, method): synchronize = True) randomise.inputs.tfce = True randomise.inputs.num_perm = 10000 - randomise.inputs. = True - randomise.inputs.pthreshold = 0.05 + randomise.inputs.c_thresh = 0.05 + group_level.connect(mask_intersection, 'out_file', randomise, 'mask') + group_level.connect(merge_copes, 'merged_file', randomise, 'in_file') + group_level.connect(specify_model, 'design_con', randomise, 'tcon') + group_level.connect(specify_model, 'design_mat', randomise, 'design_mat') - -tcon -mask -fcon -design_mat - - (estimate_model, randomise, - ('zstats', 'in_file') - ('copes', 'cope_file')) - - - - - # Datasink Node - save important files + # Datasink Node - Save important files data_sink = Node(DataSink(), name = 'data_sink') data_sink.inputs.base_directory = self.directories.output_dir - (estimate_model, data_sink, [ - ('zstats', f'group_level_analysis_{method}_nsub_{nb_subjects}.@zstats'), - ('tstats', f'group_level_analysis_{method}_nsub_{nb_subjects}.@tstats') - ]), - (cluster, data_sink, [ - ('threshold_file', f'group_level_analysis_{method}_nsub_{nb_subjects}.@thresh'), - ('pval_file', f'group_level_analysis_{method}_nsub_{nb_subjects}.@pval')]) - ]) + group_level.connect(estimate_model, 'zstats', data_sink, + f'group_level_analysis_{method}_nsub_{nb_subjects}.@zstats') + group_level.connect(estimate_model, 'tstats', data_sink, + f'group_level_analysis_{method}_nsub_{nb_subjects}.@tstats') + group_level.connect(randomise,'t_corrected_p_files', data_sink, + f'group_level_analysis_{method}_nsub_{nb_subjects}.@t_corrected_p_files') + group_level.connect(randomise,'t_p_files', data_sink, + f'group_level_analysis_{method}_nsub_{nb_subjects}.@t_p_files') if method in ('equalIndifference', 'equalRange'): # Setup a one sample t-test @@ -653,6 +643,8 @@ def get_group_level_analysis_sub_workflow(self, method): ) get_group_subjects.inputs.list_1 = get_group(method) get_group_subjects.inputs.list_2 = self.subject_list + group_level.connect(get_group_subjects, 'out_list', get_copes, 'elements') + group_level.connect(get_group_subjects, 'out_list', get_varcopes, 'elements') # Function Node get_one_sample_t_test_regressors # Get regressors in the equalRange and equalIndifference method case @@ -664,14 +656,8 @@ def get_group_level_analysis_sub_workflow(self, method): ), name = 'regressors_one_sample', ) - - # Add missing connections - group_level_analysis.connect([ - (get_group_subjects, get_copes, [('out_list', 'elements')]), - (get_group_subjects, get_varcopes, [('out_list', 'elements')]), - (get_group_subjects, regressors_one_sample, [('out_list', 'subject_list')]), - (regressors_one_sample, specify_model, [('regressors', 'regressors')]) - ]) + group_level.connect(get_group_subjects, 'out_list', regressors_one_sample, 'subject_list') + group_level.connect(regressors_one_sample, 'regressors', specify_model, 'regressors') elif method == 'groupComp': @@ -726,19 +712,15 @@ def get_group_level_analysis_sub_workflow(self, method): regressors_two_sample.inputs.subject_list = self.subject_list # Add missing connections - group_level_analysis.connect([ - (get_equal_range_subjects, regressors_two_sample, [ - ('out_list', 'equal_range_ids') - ]), - (get_equal_indifference_subjects, regressors_two_sample, [ - ('out_list', 'equal_indifference_ids') - ]), - (regressors_two_sample, specify_model, [ - ('regressors', 'regressors'), - ('groups', 'groups')]) - ]) - - return group_level_analysis + group_level.connect( + get_equal_range_subjects, 'out_list', regressors_two_sample, 'equal_range_ids') + group_level.connect( + get_equal_indifference_subjects, 'out_list', + regressors_two_sample, 'equal_indifference_ids') + group_level.connect(regressors_two_sample, 'regressors', specify_model, 'regressors') + group_level.connect(regressors_two_sample, 'groups', specify_model, 'groups') + + return group_level def get_group_level_outputs(self): """ Return all names for the files the group level analysis is supposed to generate. """ @@ -748,7 +730,7 @@ def get_group_level_outputs(self): 'contrast_id': self.contrast_list, 'method': ['equalRange', 'equalIndifference'], 'file': [ - '_cluster0/zstat1_pval.nii.gz', + '_cluster0/zstat1_pval.nii.gz', # TODO : output for randomise '_cluster0/zstat1_threshold.nii.gz', '_cluster1/zstat2_pval.nii.gz', '_cluster1/zstat2_threshold.nii.gz', @@ -772,7 +754,7 @@ def get_group_level_outputs(self): parameters = { 'contrast_id': self.contrast_list, 'file': [ - '_cluster0/zstat1_pval.nii.gz', + '_cluster0/zstat1_pval.nii.gz', # TODO : output for randomise '_cluster0/zstat1_threshold.nii.gz', 'tstat1.nii.gz', 'zstat1.nii.gz' diff --git a/tests/pipelines/test_team_B5I6.py b/tests/pipelines/test_team_B5I6.py index f0ed393f..d19aaf32 100644 --- a/tests/pipelines/test_team_B5I6.py +++ b/tests/pipelines/test_team_B5I6.py @@ -52,11 +52,11 @@ def test_outputs(): # 1 - 1 subject outputs pipeline.subject_list = ['001'] - helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*1, 4*4*1 + 4*1, 8*4*2 + 4*4, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*1, 4*6*4*1 + 4*1, 8*4*2 + 4*4, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*4, 4*4*4 + 4*4, 8*4*2 + 4*4, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*4, 4*6*4*4*4 + 4*4, 8*4*2 + 4*4, 18]) @staticmethod @mark.unit_test @@ -138,9 +138,6 @@ def test_confounds_file_no_outliers(temporary_data_dir): temporary_data_dir, confounds_node.name, 'confounds_file_sub-sid_run-rid.tsv')) assert exists(created_confounds_file) - from shutil import copy - copy(created_confounds_file, '/work/toot.tsv') - # Check contents assert cmp(reference_file, created_confounds_file) @@ -174,9 +171,6 @@ def test_confounds_file_outliers(temporary_data_dir): temporary_data_dir, confounds_node.name, 'confounds_file_sub-sid_run-rid.tsv')) assert exists(created_confounds_file) - from shutil import copy - copy(created_confounds_file, '/work/toot_outl.tsv') - # Check contents assert cmp(reference_file, created_confounds_file) From 35ab05f200eb8ad122b57e97f846da9f78f65dc4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Fri, 12 Apr 2024 11:50:41 +0200 Subject: [PATCH 05/12] Removing DC61 code and data --- narps_open/pipelines/__init__.py | 2 +- narps_open/pipelines/team_B5I6.py | 10 - narps_open/pipelines/team_DC61.py | 735 ------------------ tests/pipelines/test_team_DC61.py | 163 ---- .../pipelines/team_DC61/confounds.tsv | 3 - 5 files changed, 1 insertion(+), 912 deletions(-) delete mode 100644 narps_open/pipelines/team_DC61.py delete mode 100644 tests/pipelines/test_team_DC61.py delete mode 100644 tests/test_data/pipelines/team_DC61/confounds.tsv diff --git a/narps_open/pipelines/__init__.py b/narps_open/pipelines/__init__.py index e02d2f22..d0fc3dff 100644 --- a/narps_open/pipelines/__init__.py +++ b/narps_open/pipelines/__init__.py @@ -44,7 +44,7 @@ 'B5I6': 'PipelineTeamB5I6', 'C22U': None, 'C88N': 'PipelineTeamC88N', - 'DC61': 'PipelineTeamDC61', + 'DC61': None, 'E3B6': None, 'E6R3': None, 'I07H': None, diff --git a/narps_open/pipelines/team_B5I6.py b/narps_open/pipelines/team_B5I6.py index bdc65f3a..19539573 100644 --- a/narps_open/pipelines/team_B5I6.py +++ b/narps_open/pipelines/team_B5I6.py @@ -158,26 +158,16 @@ def get_confounds_file(filepath, subject_id, run_id): # Compute outliers of `non-stdDVARS` dvars_outliers = array(data_frame['non-stdDVARS']) - print(dvars_outliers) - percentile_75 = nanpercentile(dvars_outliers, 75) interquartile_range = percentile_75 - nanpercentile(dvars_outliers, 25) dvars_outliers = (dvars_outliers > (percentile_75 + 1.5 * interquartile_range)).astype(int) - print(dvars_outliers) - print(percentile_75 + 1.5 * interquartile_range) - # Compute outliers of `FramewiseDisplacement` fd_outliers = array(data_frame['FramewiseDisplacement']) - print(fd_outliers) - percentile_75 = nanpercentile(fd_outliers, 75) interquartile_range = percentile_75 - nanpercentile(fd_outliers, 25) fd_outliers = (fd_outliers > (percentile_75 + 1.5 * interquartile_range)).astype(int) - print(fd_outliers) - print(percentile_75 + 1.5 * interquartile_range) - # Build an outlier regressor outlier_regressor = fd_outliers | dvars_outliers diff --git a/narps_open/pipelines/team_DC61.py b/narps_open/pipelines/team_DC61.py deleted file mode 100644 index 3c5dae4e..00000000 --- a/narps_open/pipelines/team_DC61.py +++ /dev/null @@ -1,735 +0,0 @@ -#!/usr/bin/python -# coding: utf-8 - -""" Write the work of NARPS' team DC61 using Nipype """ - -from os.path import join -from itertools import product - -from nipype import Workflow, Node, MapNode -from nipype.interfaces.utility import IdentityInterface, Function -from nipype.interfaces.io import SelectFiles, DataSink -from nipype.interfaces.spm import ( - Smooth, - OneSampleTTestDesign, EstimateModel, EstimateContrast, - Level1Design, TwoSampleTTestDesign, Threshold - ) -from nipype.algorithms.modelgen import SpecifySPMModel -from nipype.algorithms.misc import Gunzip - -from narps_open.pipelines import Pipeline -from narps_open.data.task import TaskInformation -from narps_open.data.participants import get_group -from narps_open.core.interfaces import InterfaceFactory -from narps_open.core.common import list_intersection, elements_in_string, clean_list -from narps_open.utils.configuration import Configuration - -class PipelineTeamDC61(Pipeline): - """ A class that defines the pipeline of team DC61. """ - - def __init__(self): - super().__init__() - self.fwhm = 5.0 - self.team_id = 'DC61' - self.contrast_list = ['0001', '0002'] - self.subject_level_contrasts = [ - ('effect_of_gain', 'T', - [f'gamble_run{r}xgain_param^1' for r in range(1, len(self.run_list) + 1)], - [0.25]*len(self.run_list)), - ('effect_of_loss', 'T', - [f'gamble_run{r}xloss_param^1' for r in range(1, len(self.run_list) + 1)], - [0.25]*len(self.run_list)) - ] - - def get_preprocessing(self): - """ No preprocessing has been done by team DC61 """ - return None - - def get_run_level_analysis(self): - """ No run level analysis has been done by team DC61 """ - return None - - # @staticmethod # Starting python 3.10, staticmethod should be used here - # Otherwise it produces a TypeError: 'staticmethod' object is not callable - def get_subject_information(event_files: list): - """ Create Bunchs for SpecifySPMModel. - - Parameters : - - event_files: list of str, list of events files (one per run) for the subject - - Returns : - - subject_info : list of Bunch for 1st level analysis. - """ - from nipype.interfaces.base import Bunch - - subject_info = [] - - for run_id, event_file in enumerate(event_files): - - onsets = [] - durations = [] - gain_value = [] - loss_value = [] - reaction_time = [] - - # Parse the events file - with open(event_file, 'rt') as file: - next(file) # skip the header - - for line in file: - info = line.strip().split() - - onsets.append(float(info[0])) - durations.append(float(info[1])) - gain_value.append(float(info[2])) - loss_value.append(float(info[3])) - reaction_time.append(float(info[4])) - - # Create a Bunch for the run - subject_info.append( - Bunch( - conditions = [f'gamble_run{run_id + 1}'], - onsets = [onsets], - durations = [durations], - amplitudes = None, - tmod = None, - pmod = [ - Bunch( - name = ['gain_param', 'loss_param', 'rt_param'], - poly = [1, 1, 1], - param = [gain_value, loss_value, reaction_time] - ) - ], - regressor_names = None, - regressors = None - )) - - return subject_info - - # @staticmethod # Starting python 3.10, staticmethod should be used here - # Otherwise it produces a TypeError: 'staticmethod' object is not callable - def get_confounds_file(filepath, subject_id, run_id, working_dir): - """ - Create a new tsv files with only desired confounds per subject per run. - Also computes the first derivative of the motion parameters. - - Parameters : - - filepath : path to the subject confounds file - - subject_id : related subject id - - run_id : related run id - - working_dir: str, name of the directory for intermediate results - - Return : - - confounds_file : path to new file containing only desired confounds - """ - from os import makedirs - from os.path import join - - from pandas import DataFrame, read_csv - from numpy import array, transpose, diff, square, insert - - # Open original confounds file - data_frame = read_csv(filepath, sep = '\t', header=0) - - # Extract confounds we want to use for the model - retained_parameters = DataFrame(transpose(array([ - data_frame['X'], data_frame['Y'], data_frame['Z'], - data_frame['RotX'], data_frame['RotY'], data_frame['RotZ'], - insert(diff(data_frame['X']), 0, 0), - insert(diff(data_frame['Y']), 0, 0), - insert(diff(data_frame['Z']), 0, 0), - insert(diff(data_frame['RotX']), 0, 0), - insert(diff(data_frame['RotY']), 0, 0), - insert(diff(data_frame['RotZ']), 0, 0), - square(data_frame['X']), square(data_frame['Y']), square(data_frame['Z']), - square(data_frame['RotX']), square(data_frame['RotY']), square(data_frame['RotZ']), - insert(square(diff(data_frame['X'])), 0, 0), - insert(square(diff(data_frame['Y'])), 0, 0), - insert(square(diff(data_frame['Z'])), 0, 0), - insert(square(diff(data_frame['RotX'])), 0, 0), - insert(square(diff(data_frame['RotY'])), 0, 0), - insert(square(diff(data_frame['RotZ'])), 0, 0) - ]))) - - # Write confounds to a file - confounds_file = join(working_dir, 'confounds_files', - f'confounds_file_sub-{subject_id}_run-{run_id}.tsv') - - makedirs(join(working_dir, 'confounds_files'), exist_ok = True) - - with open(confounds_file, 'w', encoding = 'utf-8') as writer: - writer.write(retained_parameters.to_csv( - sep = '\t', index = False, header = False, na_rep = '0.0')) - - return confounds_file - - def get_subject_level_analysis(self): - """ - Create the subject level analysis workflow. - - Returns: - - subject_level : nipype.WorkFlow - """ - # Initialize preprocessing workflow to connect nodes along the way - subject_level = Workflow( - base_dir = self.directories.working_dir, name = 'subject_level' - ) - - # Identity interface Node - to iterate over subject_id and run - info_source = Node( - IdentityInterface(fields = ['subject_id']), - name = 'info_source') - info_source.iterables = [('subject_id', self.subject_list)] - - # Select files from derivatives - templates = { - 'func': join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', - 'sub-{subject_id}_task-MGT_run-*_bold_space-MNI152NLin2009cAsym_preproc.nii.gz'), - 'confounds' : join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', - 'sub-{subject_id}_task-MGT_run-*_bold_confounds.tsv'), - 'events': join('sub-{subject_id}', 'func', - 'sub-{subject_id}_task-MGT_run-*_events.tsv') - } - select_files = Node(SelectFiles(templates), name = 'select_files') - select_files.inputs.base_directory = self.directories.dataset_dir - select_files.inputs.sort_filelist = True - subject_level.connect(info_source, 'subject_id', select_files, 'subject_id') - - # Gunzip - gunzip files because SPM do not use .nii.gz files - gunzip = MapNode(Gunzip(), name = 'gunzip', iterfield=['in_file']) - subject_level.connect(select_files, 'func', gunzip, 'in_file') - - # Smoothing - smooth the func data - smooth = Node(Smooth(), name = 'smooth') - smooth.inputs.fwhm = self.fwhm - smooth.overwrite = False - subject_level.connect(gunzip, 'out_file', smooth, 'in_files') - - # Function node get_subject_info - get subject specific condition information - subject_info = Node(Function( - function = self.get_subject_information, - input_names = ['event_files'], - output_names = ['subject_info'] - ), - name = 'subject_info') - subject_level.connect(select_files, 'events', subject_info, 'event_files') - - # Function node get_confounds_file - Generate confounds files - confounds = MapNode(Function( - function = self.get_confounds_file, - input_names = ['filepath', 'subject_id', 'run_id', 'working_dir'], - output_names = ['confounds_file']), - name = 'confounds', iterfield = ['filepath', 'run_id']) - confounds.inputs.working_dir = self.directories.working_dir - confounds.inputs.run_id = self.run_list - subject_level.connect(info_source, 'subject_id', confounds, 'subject_id') - subject_level.connect(select_files, 'confounds', confounds, 'filepath') - - specify_model = Node(SpecifySPMModel(), name = 'specify_model') - specify_model.inputs.concatenate_runs = False - specify_model.inputs.input_units = 'secs' - specify_model.inputs.output_units = 'secs' - specify_model.inputs.time_repetition = TaskInformation()['RepetitionTime'] - specify_model.inputs.high_pass_filter_cutoff = 128 - specify_model.overwrite = False - subject_level.connect(subject_info, 'subject_info', specify_model, 'subject_info') - subject_level.connect(confounds, 'confounds_file', specify_model, 'realignment_parameters') - subject_level.connect(smooth, 'smoothed_files', specify_model, 'functional_runs') - - model_design = Node(Level1Design(), name = 'model_design') - model_design.inputs.bases = {'hrf': {'derivs': [1, 0]}} # Temporal derivatives - model_design.inputs.timing_units = 'secs' - model_design.inputs.interscan_interval = TaskInformation()['RepetitionTime'] - model_design.overwrite = False - subject_level.connect(specify_model, 'session_info', model_design, 'session_info') - - model_estimate = Node(EstimateModel(), name = 'model_estimate') - model_estimate.inputs.estimation_method = {'Classical': 1} - model_estimate.overwrite = False - subject_level.connect(model_design, 'spm_mat_file', model_estimate, 'spm_mat_file') - - contrast_estimate = Node(EstimateContrast(), name = 'contraste_estimate') - contrast_estimate.inputs.contrasts = self.subject_level_contrasts - contrast_estimate.config = {'execution': {'remove_unnecessary_outputs': False}} - contrast_estimate.overwrite = False - subject_level.connect(model_estimate, 'spm_mat_file', contrast_estimate, 'spm_mat_file') - subject_level.connect(model_estimate, 'beta_images', contrast_estimate, 'beta_images') - subject_level.connect( - model_estimate, 'residual_image', contrast_estimate, 'residual_image') - - # DataSink - store the wanted results in the wanted repository - data_sink = Node(DataSink(), name = 'data_sink') - data_sink.inputs.base_directory = self.directories.output_dir - subject_level.connect( - contrast_estimate, 'con_images', data_sink, f'{subject_level.name}.@con_images') - subject_level.connect( - contrast_estimate, 'spm_mat_file', data_sink, f'{subject_level.name}.@spm_mat_file') - - # Remove large files, if requested - if Configuration()['pipelines']['remove_unused_data']: - - # Remove Node - Remove gunzip files once they are no longer needed - remove_gunzip = MapNode( - InterfaceFactory.create('remove_parent_directory'), - name = 'remove_gunzip', - iterfield = ['file_name'] - ) - - # Remove Node - Remove smoothed files once they are no longer needed - remove_smooth = MapNode( - InterfaceFactory.create('remove_parent_directory'), - name = 'remove_smooth', - iterfield = ['file_name'] - ) - - # Add connections - subject_level.connect([ - (smooth, remove_gunzip, [('smoothed_files', '_')]), - (gunzip, remove_gunzip, [('out_file', 'file_name')]), - (data_sink, remove_smooth, [('out_file', '_')]), - (smooth, remove_smooth, [('smoothed_files', 'file_name')]) - ]) - - return subject_level - - def get_subject_level_outputs(self): - """ Return the names of the files the subject level analysis is supposed to generate. """ - - templates = [join( - self.directories.output_dir, - 'subject_level', '_subject_id_{subject_id}', f'con_{contrast_id}.nii')\ - for contrast_id in self.contrast_list] - templates += [join( - self.directories.output_dir, - 'subject_level', '_subject_id_{subject_id}', 'SPM.mat')] - templates += [join( - self.directories.output_dir, - 'subject_level', '_subject_id_{subject_id}', f'spmT_{contrast_id}.nii')\ - for contrast_id in self.contrast_list] - - # Format with subject_ids - return_list = [] - for template in templates: - return_list += [template.format(subject_id = s) for s in self.subject_list] - - return return_list - - def get_group_level_contrasts(subject_level_contrast: str): - """ Return a contrast list for EstimateContrast - - Parameters : - - subject_level_contrast: str, id of the subject level contrast : - either 'effect_of_gain' - or 'effect_of_loss' - - Returns : - - list of contrasts for the EstimateContrast interface. - """ - if subject_level_contrast == 'effect_of_gain': - return [ - ['gain_param_range', 'T', ['equalIndifference', 'equalRange'], [0, 1]], - ['gain_param_indiff', 'T', ['equalIndifference', 'equalRange'], [1, 0]] - ] - - if subject_level_contrast == 'effect_of_loss': - return [ - ['loss_param_range', 'F', ['equalIndifference', 'equalRange'], [0, 1]], - ['loss_param_indiff', 'F', ['equalIndifference', 'equalRange'], [1, 0]] - ] - - return [] - - def get_group_covariates(subjects: list): - """ Return a covariates list for OneSampleTTestDesign - - Parameters : - - subjects: list, list of the subjects involved in the design - - Returns : - - list of covariates for the OneSampleTTestDesign interface. - """ - from narps_open.data.participants import get_group - - return [ - dict( - vector = [1 if s in get_group('equalRange') else 0 for s in subjects], - name = 'equalRange'), - dict( - vector = [1 if s in get_group('equalIndifference') else 0 for s in subjects], - name = 'equalIndifference') - ] - - def get_group_level_analysis(self): - """ - Return all workflows for the group level analysis. - - Returns; - - a list of nipype.WorkFlow - """ - return [ - self.get_group_level_analysis_single_group(), - self.get_group_level_analysis_group_comparison() - ] - - def get_group_level_analysis_single_group(self): - """ - Return a workflow for the group level analysis. - - Returns: - - group_level: nipype.WorkFlow - """ - # Compute the number of participants used to do the analysis - nb_subjects = len(self.subject_list) - - # Create the group level workflow - group_level = Workflow( - base_dir = self.directories.working_dir, - name = f'group_level_analysis_nsub_{nb_subjects}') - - # IDENTITY INTERFACE - Iterate over the list of subject-level contrasts - info_source = Node(IdentityInterface(fields=['contrast_id', 'contrast_name']), - name = 'info_source') - info_source.iterables = [ - ('contrast_id', self.contrast_list), - ('contrast_name', [c[0] for c in self.subject_level_contrasts])] - info_source.synchronize = True - - # SELECT FILES - Get files from subject-level analysis - templates = { - 'contrasts': join(self.directories.output_dir, - 'subject_level', '_subject_id_*', 'con_{contrast_id}.nii') - } - select_files = Node(SelectFiles(templates), name = 'select_files') - select_files.inputs.sort_filelist = True - select_files.inputs.base_directory = self.directories.dataset_dir - - # Create a function to complete the subject ids out from the get_equal_*_subjects nodes - # If not complete, subject id '001' in search patterns - # would match all contrast files with 'con_0001.nii'. - complete_subject_ids = lambda l : [f'_subject_id_{a}' for a in l] - - # GET CONTRAST FILES - # Function Node elements_in_string - # Get contrast files for required subjects - # Note : using a MapNode with elements_in_string requires using clean_list to remove - # None values from the out_list - get_contrast_files = MapNode(Function( - function = elements_in_string, - input_names = ['input_str', 'elements'], - output_names = ['out_list'] - ), - name = 'get_contrast_files', iterfield = 'input_str' - ) - get_contrast_files.inputs.elements = complete_subject_ids(self.subject_list) - group_level.connect(select_files, 'contrasts', get_contrast_files, 'input_str') - - # GET COVARIATES - # Function Node get_group_covariates - # Get groups as covariates input for OneSampleTTestDesign - get_covariates = Node(Function( - function = self.get_group_covariates, - input_names = ['subjects'], - output_names = ['covariates'] - ), - name = 'get_covariates' - ) - get_covariates.inputs.subjects = self.subject_list - - # ONE SAMPLE T-TEST DESIGN - Create a t test design for hypothesis testing inside a group - one_sample_t_test = Node(OneSampleTTestDesign(), name = 'one_sample_t_test') - group_level.connect(get_covariates, 'covariates', one_sample_t_test, 'covariates') - group_level.connect( - get_contrast_files, ('out_list', clean_list), one_sample_t_test, 'in_files') - - # ESTIMATE MODEL - estimate the parameters of the model - # Even for second level it should be 'Classical': 1. - model_estimate = Node(EstimateModel(), name = 'model_estimate') - model_estimate.inputs.estimation_method = {'Classical': 1} - group_level.connect(one_sample_t_test, 'spm_mat_file', model_estimate, 'spm_mat_file') - - # Function Node get_group_level_contrasts - # Get a contrast list for EstimateContrast - get_contrasts = Node(Function( - function = self.get_group_level_contrasts, - input_names = ['subject_level_contrast'], - output_names = ['contrasts'] - ), - name = 'get_contrasts' - ) - - # ESTIMATE CONTRASTS - estimates simple group contrast - contrast_estimate = Node(EstimateContrast(), name = 'contrast_estimate') - contrast_estimate.inputs.group_contrast = True - group_level.connect(get_contrasts, 'contrasts', contrast_estimate, 'contrasts') - group_level.connect(model_estimate, 'spm_mat_file', contrast_estimate, 'spm_mat_file') - group_level.connect(model_estimate, 'beta_images', contrast_estimate, 'beta_images') - group_level.connect(model_estimate, 'residual_image', contrast_estimate, 'residual_image') - - # THRESHOLD - Create thresholded maps - threshold = MapNode(Threshold(), name = 'threshold', - iterfield = ['stat_image', 'contrast_index']) - threshold.inputs.use_fwe_correction = True - threshold.inputs.height_threshold_type = 'p-value' - threshold.inputs.force_activation = False - threshold.inputs.height_threshold = 0.05 - threshold.inputs.contrast_index = [1, 2] - group_level.connect(contrast_estimate, 'spm_mat_file', threshold, 'spm_mat_file') - group_level.connect(contrast_estimate, 'spmT_images', threshold, 'stat_image') - - # Datasink - save important files - data_sink = Node(DataSink(), name = 'data_sink') - data_sink.inputs.base_directory = self.directories.output_dir - group_level.connect([ - (model_estimate, data_sink, [ - ('mask_image', f'{group_level.name}.@mask')]), - (contrast_estimate, data_sink, [ - ('spm_mat_file', f'{group_level.name}.@spm_mat'), - ('spmT_images', f'{group_level.name}.@T'), - ('con_images', f'{group_level.name}.@con')]), - (threshold, data_sink, [ - ('thresholded_map', f'{group_level.name}.@thresh')]) - ]) - - return group_level - - def get_group_level_analysis_group_comparison(self): - """ - Return a workflow for the group level analysis in the group comparison case. - - Returns: - - group_level_analysis: nipype.WorkFlow - """ - # Compute the number of participants used to do the analysis - nb_subjects = len(self.subject_list) - - # Create the group level workflow - group_level_analysis = Workflow( - base_dir = self.directories.working_dir, - name = f'group_level_analysis_groupComp_nsub_{nb_subjects}') - - # SELEC FILES - Get files from subject level analysis (effect_of_loss contrast only) - templates = { - 'contrasts': join(self.directories.output_dir, - 'subject_level', '_subject_id_*', 'con_0002.nii') - } - select_files = Node(SelectFiles(templates), name = 'select_files') - select_files.inputs.sort_filelist = True - select_files.inputs.base_directory = self.directories.dataset_dir - - # GET SUBJECT LIST IN EACH GROUP - # Function Node get_group_subjects - # Get subjects in the group and in the subject_list - get_equal_indifference_subjects = Node(Function( - function = list_intersection, - input_names = ['list_1', 'list_2'], - output_names = ['out_list'] - ), - name = 'get_equal_indifference_subjects' - ) - get_equal_indifference_subjects.inputs.list_1 = get_group('equalIndifference') - get_equal_indifference_subjects.inputs.list_2 = self.subject_list - - # Function Node get_group_subjects - # Get subjects in the group and in the subject_list - get_equal_range_subjects = Node(Function( - function = list_intersection, - input_names = ['list_1', 'list_2'], - output_names = ['out_list'] - ), - name = 'get_equal_range_subjects' - ) - get_equal_range_subjects.inputs.list_1 = get_group('equalRange') - get_equal_range_subjects.inputs.list_2 = self.subject_list - - # Create a function to complete the subject ids out from the get_equal_*_subjects nodes - # If not complete, subject id '001' in search patterns - # would match all contrast files with 'con_0001.nii'. - complete_subject_ids = lambda l : [f'_subject_id_{a}' for a in l] - - # GET CONTRAST FILES - # Function Node elements_in_string - # Get contrast files for required subjects - # Note : using a MapNode with elements_in_string requires using clean_list to remove - # None values from the out_list - get_equal_indifference_contrasts = MapNode(Function( - function = elements_in_string, - input_names = ['input_str', 'elements'], - output_names = ['out_list'] - ), - name = 'get_equal_indifference_contrasts', iterfield = 'input_str' - ) - group_level_analysis.connect( - select_files, 'contrasts', get_equal_indifference_contrasts, 'input_str') - group_level_analysis.connect( - get_equal_indifference_subjects, ('out_list', complete_subject_ids), - get_equal_indifference_contrasts, 'elements') - - get_equal_range_contrasts = MapNode(Function( - function = elements_in_string, - input_names = ['input_str', 'elements'], - output_names = ['out_list'] - ), - name = 'get_equal_range_contrasts', iterfield = 'input_str' - ) - group_level_analysis.connect( - select_files, 'contrasts', get_equal_range_contrasts, 'input_str') - group_level_analysis.connect( - get_equal_range_subjects, ('out_list', complete_subject_ids), - get_equal_range_contrasts, 'elements') - - # TWO SAMPLE T-TEST DESIGN - Create a t test design for group comparison - two_sample_t_test = Node(TwoSampleTTestDesign(), name = 'two_sample_t_test') - group_level_analysis.connect( - get_equal_range_contrasts, ('out_list', clean_list), - two_sample_t_test, 'group1_files') - group_level_analysis.connect( - get_equal_indifference_contrasts, ('out_list', clean_list), - two_sample_t_test, 'group2_files') - - # ESTIMATE MODEL - Estimate the parameters of the model - # Even for second level it should be 'Classical': 1. - model_estimate = Node(EstimateModel(), name = 'model_estimate') - model_estimate.inputs.estimation_method = {'Classical': 1} - group_level_analysis.connect( - two_sample_t_test, 'spm_mat_file', model_estimate, 'spm_mat_file') - - # ESTIMATE CONTRASTS - Estimates contrasts - contrast_estimate = Node(EstimateContrast(), name = 'contrast_estimate') - contrast_estimate.inputs.group_contrast = True - contrast_estimate.inputs.contrasts = [ - ['Eq range vs Eq indiff in loss', 'T', ['Group_{1}', 'Group_{2}'], [1, -1]] - ] - group_level_analysis.connect([ - (model_estimate, contrast_estimate, [ - ('spm_mat_file', 'spm_mat_file'), - ('beta_images', 'beta_images'), - ('residual_image', 'residual_image') - ])]) - - # THRESHOLD - Create thresholded maps - threshold = Node(Threshold(), name = 'threshold') - threshold.inputs.use_fwe_correction = True - threshold.inputs.height_threshold_type = 'p-value' - threshold.inputs.force_activation = False - threshold.inputs.height_threshold = 0.05 - threshold.inputs.contrast_index = 1 - group_level_analysis.connect([ - (contrast_estimate, threshold, [ - ('spm_mat_file', 'spm_mat_file'), - ('spmT_images', 'stat_image') - ])]) - - # DATA SINK - save important files - data_sink = Node(DataSink(), name = 'data_sink') - data_sink.inputs.base_directory = self.directories.output_dir - group_level_analysis.connect([ - (model_estimate, data_sink, [ - ('mask_image', f'{group_level_analysis.name}.@mask')]), - (contrast_estimate, data_sink, [ - ('spm_mat_file', f'{group_level_analysis.name}.@spm_mat'), - ('spmT_images', f'{group_level_analysis.name}.@T'), - ('con_images', f'{group_level_analysis.name}.@con')]), - (threshold, data_sink, [ - ('thresholded_map', f'{group_level_analysis.name}.@thresh')]) - ]) - - return group_level_analysis - - def get_group_level_outputs(self): - """ Return all names for the files the group level analysis is supposed to generate. """ - - # Handle equalRange and equalIndifference - parameters = { - 'contrast_id': self.contrast_list, - 'method': ['equalRange', 'equalIndifference'], - 'file': [ - 'con_0001.nii', 'con_0002.nii', 'mask.nii', 'SPM.mat', - 'spmT_0001.nii', 'spmT_0002.nii', - join('_threshold0', 'spmT_0001_thr.nii'), join('_threshold1', 'spmT_0002_thr.nii') - ], - 'nb_subjects' : [str(len(self.subject_list))] - } - - parameter_sets = product(*parameters.values()) - template = join( - self.directories.output_dir, - 'group_level_analysis_{method}_nsub_{nb_subjects}', - '_contrast_id_{contrast_id}', - '{file}' - ) - return_list = [template.format(**dict(zip(parameters.keys(), parameter_values)))\ - for parameter_values in parameter_sets] - - # Handle groupComp - parameters = { - 'contrast_id': self.contrast_list, - 'method': ['groupComp'], - 'file': [ - 'con_0001.nii', 'mask.nii', 'SPM.mat', 'spmT_0001.nii', - join('_threshold0', 'spmT_0001_thr.nii') - ], - 'nb_subjects' : [str(len(self.subject_list))] - } - parameter_sets = product(*parameters.values()) - template = join( - self.directories.output_dir, - 'group_level_analysis_{method}_nsub_{nb_subjects}', - '_contrast_id_{contrast_id}', - '{file}' - ) - return_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ - for parameter_values in parameter_sets] - - return return_list - - def get_hypotheses_outputs(self): - """ Return all hypotheses output file names. """ - nb_sub = len(self.subject_list) - files = [ - # Hypothesis 1 - join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', - '_contrast_id_0001', '_threshold0', 'spmT_0001_thr.nii'), - join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', - '_contrast_id_0001', 'spmT_0001.nii'), - # Hypothesis 2 - join(f'group_level_analysis_equalRange_nsub_{nb_sub}', - '_contrast_id_0001', '_threshold0', 'spmT_0001_thr.nii'), - join(f'group_level_analysis_equalRange_nsub_{nb_sub}', - '_contrast_id_0001', 'spmT_0001.nii'), - # Hypothesis 3 - join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', - '_contrast_id_0001', '_threshold0', 'spmT_0001_thr.nii'), - join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', - '_contrast_id_0001', 'spmT_0001.nii'), - # Hypothesis 4 - join(f'group_level_analysis_equalRange_nsub_{nb_sub}', - '_contrast_id_0001', '_threshold0', 'spmT_0001_thr.nii'), - join(f'group_level_analysis_equalRange_nsub_{nb_sub}', - '_contrast_id_0001', 'spmT_0001.nii'), - # Hypothesis 5 - join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', - '_contrast_id_0002', '_threshold1', 'spmT_0002_thr.nii'), - join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', - '_contrast_id_0002', 'spmT_0002.nii'), - # Hypothesis 6 - join(f'group_level_analysis_equalRange_nsub_{nb_sub}', - '_contrast_id_0002', '_threshold1', 'spmT_0002_thr.nii'), - join(f'group_level_analysis_equalRange_nsub_{nb_sub}', - '_contrast_id_0002', 'spmT_0002.nii'), - # Hypothesis 7 - join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', - '_contrast_id_0002', '_threshold0', 'spmT_0001_thr.nii'), - join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', - '_contrast_id_0002', 'spmT_0001.nii'), - # Hypothesis 8 - join(f'group_level_analysis_equalRange_nsub_{nb_sub}', - '_contrast_id_0002', '_threshold0', 'spmT_0001_thr.nii'), - join(f'group_level_analysis_equalRange_nsub_{nb_sub}', - '_contrast_id_0002', 'spmT_0001.nii'), - # Hypothesis 9 - join(f'group_level_analysis_groupComp_nsub_{nb_sub}', - '_contrast_id_0002', '_threshold0', 'spmT_0001_thr.nii'), - join(f'group_level_analysis_groupComp_nsub_{nb_sub}', - '_contrast_id_0002', 'spmT_0001.nii') - ] - return [join(self.directories.output_dir, f) for f in files] diff --git a/tests/pipelines/test_team_DC61.py b/tests/pipelines/test_team_DC61.py deleted file mode 100644 index d25a81e3..00000000 --- a/tests/pipelines/test_team_DC61.py +++ /dev/null @@ -1,163 +0,0 @@ -#!/usr/bin/python -# coding: utf-8 - -""" Tests of the 'narps_open.pipelines.team_DC61' module. - -Launch this test with PyTest - -Usage: -====== - pytest -q test_team_DC61.py - pytest -q test_team_DC61.py -k -""" -from os.path import join, exists -from filecmp import cmp - -from pytest import helpers, mark, fixture -from nipype import Workflow -from nipype.interfaces.base import Bunch - -from narps_open.utils.configuration import Configuration -from narps_open.pipelines.team_DC61 import PipelineTeamDC61 - -@fixture -def mock_participants_data(mocker): - """ A fixture to provide mocked data from the test_data directory """ - - mocker.patch( - 'narps_open.data.participants.Configuration', - return_value = { - 'directories': { - 'dataset': join( - Configuration()['directories']['test_data'], - 'data', 'participants') - } - } - ) - -class TestPipelinesTeamDC61: - """ A class that contains all the unit tests for the PipelineTeamDC61 class.""" - - @staticmethod - @mark.unit_test - def test_create(): - """ Test the creation of a PipelineTeamDC61 object """ - - pipeline = PipelineTeamDC61() - - # 1 - check the parameters - assert pipeline.fwhm == 5.0 - assert pipeline.team_id == 'DC61' - - # 2 - check workflows - assert pipeline.get_preprocessing() is None - assert pipeline.get_run_level_analysis() is None - assert isinstance(pipeline.get_subject_level_analysis(), Workflow) - group_level = pipeline.get_group_level_analysis() - assert len(group_level) == 2 - for sub_workflow in group_level: - assert isinstance(sub_workflow, Workflow) - - @staticmethod - @mark.unit_test - def test_outputs(): - """ Test the expected outputs of a PipelineTeamDC61 object """ - pipeline = PipelineTeamDC61() - # 1 - 1 subject outputs - pipeline.subject_list = ['001'] - helpers.test_pipeline_outputs(pipeline, [0, 0, 5, 8*2*2 + 5*2, 18]) - - # 2 - 4 subjects outputs - pipeline.subject_list = ['001', '002', '003', '004'] - helpers.test_pipeline_outputs(pipeline, [0, 0, 20, 8*2*2 + 5*2, 18]) - - @staticmethod - @mark.unit_test - def test_subject_information(): - """ Test the get_subject_information method """ - - # Get test files - test_file = join(Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') - info = PipelineTeamDC61.get_subject_information([test_file, test_file]) - - # Compare bunches to expected - bunch = info[0] - assert isinstance(bunch, Bunch) - assert bunch.conditions == ['gamble_run1'] - helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 19.535, 27.535, 36.435]]) - helpers.compare_float_2d_arrays(bunch.durations, [[4.0, 4.0, 4.0, 4.0, 4.0]]) - assert bunch.amplitudes is None - assert bunch.tmod is None - assert bunch.pmod[0].name == ['gain_param', 'loss_param', 'rt_param'] - assert bunch.pmod[0].poly == [1, 1, 1] - helpers.compare_float_2d_arrays(bunch.pmod[0].param, [ - [14.0, 34.0, 38.0, 10.0, 16.0], - [6.0, 14.0, 19.0, 15.0, 17.0], - [2.388, 2.289, 0.0, 2.08, 2.288] - ]) - assert bunch.regressor_names is None - assert bunch.regressors is None - - bunch = info[1] - assert isinstance(bunch, Bunch) - assert bunch.conditions == ['gamble_run2'] - helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 19.535, 27.535, 36.435]]) - helpers.compare_float_2d_arrays(bunch.durations, [[4.0, 4.0, 4.0, 4.0, 4.0]]) - assert bunch.amplitudes is None - assert bunch.tmod is None - assert bunch.pmod[0].name == ['gain_param', 'loss_param', 'rt_param'] - assert bunch.pmod[0].poly == [1, 1, 1] - helpers.compare_float_2d_arrays(bunch.pmod[0].param, [ - [14.0, 34.0, 38.0, 10.0, 16.0], - [6.0, 14.0, 19.0, 15.0, 17.0], - [2.388, 2.289, 0.0, 2.08, 2.288] - ]) - assert bunch.regressor_names is None - assert bunch.regressors is None - - @staticmethod - @mark.unit_test - def test_confounds_file(temporary_data_dir): - """ Test the get_confounds_file method """ - - confounds_file = join( - Configuration()['directories']['test_data'], 'pipelines', 'confounds.tsv') - reference_file = join( - Configuration()['directories']['test_data'], 'pipelines', 'team_DC61', 'confounds.tsv') - - # Get new confounds file - PipelineTeamDC61.get_confounds_file(confounds_file, 'sid', 'rid', temporary_data_dir) - - # Check confounds file was created - created_confounds_file = join( - temporary_data_dir, 'confounds_files', 'confounds_file_sub-sid_run-rid.tsv') - assert exists(created_confounds_file) - - # Check contents - assert cmp(reference_file, created_confounds_file) - - @staticmethod - @mark.unit_test - def test_group_level_contrasts(): - """ Test the get_group_contrasts method """ - - @staticmethod - @mark.unit_test - def test_group_covariates(mock_participants_data): - """ Test the get_group_covariates method """ - subjects = ['001', '002', '003', '004'] - assert PipelineTeamDC61.get_group_covariates(subjects) == [ - dict(vector = [0, 1, 0, 1], name = 'equalRange'), - dict(vector = [1, 0, 1, 0], name = 'equalIndifference') - ] - subjects = ['001', '003', '002', '004'] - assert PipelineTeamDC61.get_group_covariates(subjects) == [ - dict(vector = [0, 0, 1, 1], name = 'equalRange'), - dict(vector = [1, 1, 0, 0], name = 'equalIndifference') - ] - - @staticmethod - @mark.pipeline_test - def test_execution(): - """ Test the execution of a PipelineTeamDC61 and compare results """ - helpers.test_pipeline_evaluation('DC61') diff --git a/tests/test_data/pipelines/team_DC61/confounds.tsv b/tests/test_data/pipelines/team_DC61/confounds.tsv deleted file mode 100644 index 1453081d..00000000 --- a/tests/test_data/pipelines/team_DC61/confounds.tsv +++ /dev/null @@ -1,3 +0,0 @@ -0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 --0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 -0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 9.937996410250001e-05 0.00098247141136 9.0559466761e-12 1.7605839969e-06 1.4760426124900002e-07 2.82878761e-08 9.937996410250001e-05 0.00098247141136 9.0559466761e-12 1.7605839969e-06 1.4760426124900002e-07 2.82878761e-08 --2.56954e-05 -0.00923735 0.0549667 0.000997278 -0.00019745 -0.000398988 0.009943254600000001 0.022107050000000003 0.05496970931 -0.00032959199999999986 0.000186743 -0.00023079800000000002 6.6025358116e-10 8.53286350225e-05 0.00302133810889 9.94563409284e-07 3.89865025e-08 1.59191424144e-07 9.886831204042119e-05 0.0004887216597025001 0.003021668941625901 1.0863088646399991e-07 3.4872948049e-08 5.326771680400001e-08 From 816a70c284a65a818eb5da35bca5a58ef25a0741 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Fri, 12 Apr 2024 15:07:39 +0200 Subject: [PATCH 06/12] Some PEP8 --- narps_open/pipelines/team_B5I6.py | 44 ++++++++++++------------- tests/pipelines/test_team_B5I6.py | 55 +++++++++++++++++-------------- 2 files changed, 51 insertions(+), 48 deletions(-) diff --git a/narps_open/pipelines/team_B5I6.py b/narps_open/pipelines/team_B5I6.py index 19539573..afd296d3 100644 --- a/narps_open/pipelines/team_B5I6.py +++ b/narps_open/pipelines/team_B5I6.py @@ -65,8 +65,10 @@ def get_subject_information(event_file): onsets_missed = [] durations_trial = [] durations_missed = [] - weights_gain = [] - weights_loss = [] + amplitudes_trial = [] + amplitudes_missed = [] + amplitudes_gain = [] + amplitudes_loss = [] with open(event_file, 'rt') as file: next(file) # skip the header @@ -77,41 +79,35 @@ def get_subject_information(event_file): if 'NoResp' in info[5]: onsets_missed.append(float(info[0])) durations_missed.append(float(info[1])) + amplitudes_missed.append(1.0) else: onsets_trial.append(float(info[0])) durations_trial.append(float(info[1])) - weights_gain.append(float(info[2])) - weights_loss.append(float(info[3])) + amplitudes_trial.append(1.0) + amplitudes_gain.append(float(info[2])) + amplitudes_loss.append(float(info[3])) # Mean center magnitudes - weights_gain = weights_gain - mean(weights_gain) - weights_loss = weights_loss - mean(weights_loss) + amplitudes_gain = amplitudes_gain - mean(amplitudes_gain) + amplitudes_loss = amplitudes_loss - mean(amplitudes_loss) # Add missed trials if any - conditions = ['trial'] - onsets = [onsets_trial] - durations = [durations_trial] + conditions = ['trial', 'gain', 'loss'] + onsets = [onsets_trial, onsets_trial, onsets_trial] + durations = [durations_trial, durations_trial, durations_trial] + amplitudes = [amplitudes_trial, amplitudes_gain, amplitudes_loss] if onsets_missed: conditions.append('missed') onsets.append(onsets_missed) durations.append(durations_missed) + amplitudes.append(amplitudes_missed) return [ Bunch( conditions = conditions, onsets = onsets, durations = durations, - amplitudes = None, - tmod = None, - pmod = [ - Bunch( - name = ['gain', 'loss'], - poly = [1, 1], - param = [weights_gain, weights_loss] - ) - ], - regressor_names = None, - regressors = None + amplitudes = amplitudes ) ] @@ -276,8 +272,10 @@ def get_run_level_analysis(self): data_sink = Node(DataSink(), name = 'data_sink') data_sink.inputs.base_directory = self.directories.output_dir run_level.connect(model_estimate, 'results_dir', data_sink, 'run_level_analysis.@results') - run_level.connect(model_generation, 'design_file', data_sink, 'run_level_analysis.@design_file') - run_level.connect(model_generation, 'design_image', data_sink, 'run_level_analysis.@design_img') + run_level.connect( + model_generation, 'design_file', data_sink, 'run_level_analysis.@design_file') + run_level.connect( + model_generation, 'design_image', data_sink, 'run_level_analysis.@design_img') # Remove large files, if requested if Configuration()['pipelines']['remove_unused_data']: @@ -607,7 +605,7 @@ def get_group_level_analysis_sub_workflow(self, method): # Datasink Node - Save important files data_sink = Node(DataSink(), name = 'data_sink') data_sink.inputs.base_directory = self.directories.output_dir - group_level.connect(estimate_model, 'zstats', data_sink, + group_level.connect(estimate_model, 'zstats', data_sink, f'group_level_analysis_{method}_nsub_{nb_subjects}.@zstats') group_level.connect(estimate_model, 'tstats', data_sink, f'group_level_analysis_{method}_nsub_{nb_subjects}.@tstats') diff --git a/tests/pipelines/test_team_B5I6.py b/tests/pipelines/test_team_B5I6.py index d19aaf32..0895b2f5 100644 --- a/tests/pipelines/test_team_B5I6.py +++ b/tests/pipelines/test_team_B5I6.py @@ -75,36 +75,41 @@ def test_subject_information(): # Compare bunches to expected bunch = info_missed[0] assert isinstance(bunch, Bunch) - assert bunch.conditions == ['trial', 'missed'] - helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 27.535, 36.435], [19.535]]) - helpers.compare_float_2d_arrays(bunch.durations, [[4.0, 4.0, 4.0, 4.0], [4.0]]) - assert bunch.amplitudes is None - assert bunch.tmod is None - assert bunch.regressor_names is None - assert bunch.regressors is None - bunch = bunch.pmod[0] - assert isinstance(bunch, Bunch) - assert bunch.name == ['gain', 'loss'] - assert bunch.poly == [1, 1] - helpers.compare_float_2d_arrays(bunch.param, [ + assert bunch.conditions == ['trial', 'gain', 'loss', 'missed'] + helpers.compare_float_2d_arrays(bunch.onsets, [ + [4.071, 11.834, 27.535, 36.435], + [4.071, 11.834, 27.535, 36.435], + [4.071, 11.834, 27.535, 36.435], + [19.535] + ]) + helpers.compare_float_2d_arrays(bunch.durations, [ + [4.0, 4.0, 4.0, 4.0], + [4.0, 4.0, 4.0, 4.0], + [4.0, 4.0, 4.0, 4.0], + [4.0] + ]) + helpers.compare_float_2d_arrays(bunch.amplitudes, [ + [1.0, 1.0, 1.0, 1.0], [-4.5, 15.5, -8.5, -2.5], - [-7.0, 1.0, 2.0, 4.0] + [-7.0, 1.0, 2.0, 4.0], + [1.0] ]) bunch = info_no_missed[0] assert isinstance(bunch, Bunch) - assert bunch.conditions == ['trial'] - helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 27.535, 36.435]]) - helpers.compare_float_2d_arrays(bunch.durations, [[4.0, 4.0, 4.0, 4.0]]) - assert bunch.amplitudes is None - assert bunch.tmod is None - assert bunch.regressor_names is None - assert bunch.regressors is None - bunch = bunch.pmod[0] - assert isinstance(bunch, Bunch) - assert bunch.name == ['gain', 'loss'] - assert bunch.poly == [1, 1] - helpers.compare_float_2d_arrays(bunch.param, [ + assert bunch.conditions == ['trial', 'gain', 'loss'] + helpers.compare_float_2d_arrays(bunch.onsets, [ + [4.071, 11.834, 27.535, 36.435], + [4.071, 11.834, 27.535, 36.435], + [4.071, 11.834, 27.535, 36.435] + ]) + helpers.compare_float_2d_arrays(bunch.durations, [ + [4.0, 4.0, 4.0, 4.0], + [4.0, 4.0, 4.0, 4.0], + [4.0, 4.0, 4.0, 4.0] + ]) + helpers.compare_float_2d_arrays(bunch.amplitudes, [ + [1.0, 1.0, 1.0, 1.0], [-4.5, 15.5, -8.5, -2.5], [-7.0, 1.0, 2.0, 4.0] ]) From af608a648db589738a0c029ad46576fe3e9e75b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Mon, 15 Apr 2024 13:35:15 +0200 Subject: [PATCH 07/12] Randomise in group level --- narps_open/pipelines/team_B5I6.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/narps_open/pipelines/team_B5I6.py b/narps_open/pipelines/team_B5I6.py index afd296d3..324ebc24 100644 --- a/narps_open/pipelines/team_B5I6.py +++ b/narps_open/pipelines/team_B5I6.py @@ -592,7 +592,7 @@ def get_group_level_analysis_sub_workflow(self, method): # Randomise Node - Perform clustering on statistical output randomise = MapNode(Randomise(), name = 'randomise', - iterfield = ['in_file', 'dlh', 'volume', 'cope_file'], + iterfield = ['in_file', 'cope_file'], synchronize = True) randomise.inputs.tfce = True randomise.inputs.num_perm = 10000 From f054905ff35ff13ce63c5a112e2bf9f014b37b4e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Mon, 15 Apr 2024 14:13:48 +0200 Subject: [PATCH 08/12] Working on masks --- narps_open/pipelines/team_B5I6.py | 41 +++++++++++++++---------------- tests/pipelines/test_team_B5I6.py | 4 +-- 2 files changed, 22 insertions(+), 23 deletions(-) diff --git a/narps_open/pipelines/team_B5I6.py b/narps_open/pipelines/team_B5I6.py index 324ebc24..84e53a20 100644 --- a/narps_open/pipelines/team_B5I6.py +++ b/narps_open/pipelines/team_B5I6.py @@ -384,8 +384,6 @@ def get_subject_level_analysis(self): # DataSink Node - store the wanted results in the wanted directory data_sink = Node(DataSink(), name = 'data_sink') data_sink.inputs.base_directory = self.directories.output_dir - subject_level.connect( - mask_intersection, 'out_file', data_sink, 'subject_level_analysis.@mask') subject_level.connect(estimate_model, 'zstats', data_sink, 'subject_level_analysis.@stats') subject_level.connect( estimate_model, 'tstats', data_sink, 'subject_level_analysis.@tstats') @@ -411,19 +409,6 @@ def get_subject_level_outputs(self): return_list = [template.format(**dict(zip(parameters.keys(), parameter_values)))\ for parameter_values in parameter_sets] - parameters = { - 'contrast_id' : self.contrast_list, - 'subject_id' : self.subject_list, - } - parameter_sets = product(*parameters.values()) - template = join( - self.directories.output_dir, - 'subject_level_analysis', '_contrast_id_{contrast_id}_subject_id_{subject_id}', - 'sub-{subject_id}_task-MGT_run-01_bold_space-MNI152NLin2009cAsym_preproc_brain_mask_maths.nii.gz' - ) - return_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ - for parameter_values in parameter_sets] - return return_list def get_one_sample_t_test_regressors(subject_list: list) -> dict: @@ -519,9 +504,8 @@ def get_group_level_analysis_sub_workflow(self, method): 'varcope' : join(self.directories.output_dir, 'subject_level_analysis', '_contrast_id_{contrast_id}_subject_id_*', 'varcope1.nii.gz'), - 'masks': join(self.directories.output_dir, - 'subject_level_analysis', '_contrast_id_1_subject_id_*', - 'sub-*_task-MGT_run-*_bold_space-MNI152NLin2009cAsym_preproc_brain_mask_maths.nii.gz') + 'masks' : join('derivatives', 'fmriprep', 'sub-*', 'func', + 'sub-*_task-MGT_run-*_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz') } select_files = Node(SelectFiles(templates), name = 'select_files') select_files.inputs.base_directory = self.directories.results_dir @@ -553,6 +537,19 @@ def get_group_level_analysis_sub_workflow(self, method): ) group_level.connect(select_files, 'varcope', get_varcopes, 'input_str') + # Function Node elements_in_string + # Get masks for these subjects + # Note : using a MapNode with elements_in_string requires using clean_list to remove + # None values from the out_list + get_masks = MapNode(Function( + function = elements_in_string, + input_names = ['input_str', 'elements'], + output_names = ['out_list'] + ), + name = 'get_masks', iterfield = 'input_str' + ) + group_level.connect(select_files, 'masks', get_masks, 'input_str') + # Merge Node - Merge cope files merge_copes = Node(Merge(), name = 'merge_copes') merge_copes.inputs.dimension = 't' @@ -565,14 +562,14 @@ def get_group_level_analysis_sub_workflow(self, method): # Split Node - Split mask list to serve them as inputs of the MultiImageMaths node. split_masks = Node(Split(), name = 'split_masks') - split_masks.inputs.splits = [1, len(self.subject_list) - 1] + split_masks.inputs.splits = [1, len(self.subject_list * self.run_list) - 1] split_masks.inputs.squeeze = True # Unfold one-element splits removing the list - group_level.connect(select_files, 'masks', split_masks, 'inlist') + group_level.connect(get_masks, ('out_list', clean_list), split_masks, 'inlist') # MultiImageMaths Node - Create a subject mask by # computing the intersection of all run masks. mask_intersection = Node(MultiImageMaths(), name = 'mask_intersection') - mask_intersection.inputs.op_string = '-mul %s ' * (len(self.subject_list) - 1) + mask_intersection.inputs.op_string = '-add %s ' * (len(self.subject_list) - 1) + '' group_level.connect(split_masks, 'out1', mask_intersection, 'in_file') group_level.connect(split_masks, 'out2', mask_intersection, 'operand_files') @@ -633,6 +630,7 @@ def get_group_level_analysis_sub_workflow(self, method): get_group_subjects.inputs.list_2 = self.subject_list group_level.connect(get_group_subjects, 'out_list', get_copes, 'elements') group_level.connect(get_group_subjects, 'out_list', get_varcopes, 'elements') + group_level.connect(get_group_subjects, 'out_list', get_masks, 'elements') # Function Node get_one_sample_t_test_regressors # Get regressors in the equalRange and equalIndifference method case @@ -653,6 +651,7 @@ def get_group_level_analysis_sub_workflow(self, method): # Indeed the SelectFiles node asks for all (*) subjects available get_copes.inputs.elements = self.subject_list get_varcopes.inputs.elements = self.subject_list + get_masks.inputs.elements = self.subject_list # Setup a two sample t-test specify_model.inputs.contrasts = [ diff --git a/tests/pipelines/test_team_B5I6.py b/tests/pipelines/test_team_B5I6.py index 0895b2f5..956ee9ec 100644 --- a/tests/pipelines/test_team_B5I6.py +++ b/tests/pipelines/test_team_B5I6.py @@ -52,11 +52,11 @@ def test_outputs(): # 1 - 1 subject outputs pipeline.subject_list = ['001'] - helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*1, 4*6*4*1 + 4*1, 8*4*2 + 4*4, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*1, 4*6*1, 8*4*2 + 4*4, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*4, 4*6*4*4*4 + 4*4, 8*4*2 + 4*4, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*4, 4*6*4, 8*4*2 + 4*4, 18]) @staticmethod @mark.unit_test From bfaab327b83efee016c0530cd710b81afbd3744a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Mon, 15 Apr 2024 14:53:41 +0200 Subject: [PATCH 09/12] Mask computation for group level --- narps_open/pipelines/team_B5I6.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/narps_open/pipelines/team_B5I6.py b/narps_open/pipelines/team_B5I6.py index 84e53a20..4a66be47 100644 --- a/narps_open/pipelines/team_B5I6.py +++ b/narps_open/pipelines/team_B5I6.py @@ -562,14 +562,16 @@ def get_group_level_analysis_sub_workflow(self, method): # Split Node - Split mask list to serve them as inputs of the MultiImageMaths node. split_masks = Node(Split(), name = 'split_masks') - split_masks.inputs.splits = [1, len(self.subject_list * self.run_list) - 1] + split_masks.inputs.splits = [1, (len(self.subject_list) * len(self.run_list)) - 1] split_masks.inputs.squeeze = True # Unfold one-element splits removing the list group_level.connect(get_masks, ('out_list', clean_list), split_masks, 'inlist') # MultiImageMaths Node - Create a subject mask by - # computing the intersection of all run masks. + # computing the intersection of all run masks + # (all voxels included in at least 80% of masks across all runs). mask_intersection = Node(MultiImageMaths(), name = 'mask_intersection') - mask_intersection.inputs.op_string = '-add %s ' * (len(self.subject_list) - 1) + '' + mask_intersection.inputs.op_string = '-add %s ' * (len(self.subject_list) - 1)\ + + f' -thr {len(self.subject_list) * len(self.run_list) * 0.8} -bin' group_level.connect(split_masks, 'out1', mask_intersection, 'in_file') group_level.connect(split_masks, 'out2', mask_intersection, 'operand_files') From 0e7504bb68428c6bc9fdf96ee104637b0b584831 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Fri, 28 Jun 2024 09:28:25 +0200 Subject: [PATCH 10/12] Mask intersection method changed [skip ci] --- narps_open/pipelines/team_B5I6.py | 41 ++++++++++++++----------------- 1 file changed, 18 insertions(+), 23 deletions(-) diff --git a/narps_open/pipelines/team_B5I6.py b/narps_open/pipelines/team_B5I6.py index 4a66be47..2fe73a97 100644 --- a/narps_open/pipelines/team_B5I6.py +++ b/narps_open/pipelines/team_B5I6.py @@ -15,7 +15,7 @@ FSLCommand, Randomise ) from nipype.algorithms.modelgen import SpecifyModel -from nipype.interfaces.fsl.maths import MultiImageMaths +from nipype.interfaces.fsl.maths import MathsCommand from narps_open.utils.configuration import Configuration from narps_open.pipelines import Pipeline @@ -337,7 +337,7 @@ def get_subject_level_analysis(self): 'run_level_analysis', '_run_id_*_subject_id_{subject_id}', 'results', 'varcope{contrast_id}.nii.gz'), 'masks' : join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', - 'sub-{subject_id}_task-MGT_run-{run_id}_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz') + 'sub-{subject_id}_task-MGT_run-*_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz') } select_files = Node(SelectFiles(templates), name = 'select_files') select_files.inputs.base_directory= self.directories.results_dir @@ -354,18 +354,16 @@ def get_subject_level_analysis(self): merge_varcopes.inputs.dimension = 't' subject_level.connect(select_files, 'varcope', merge_varcopes, 'in_files') - # Split Node - Split mask list to serve them as inputs of the MultiImageMaths node. - split_masks = Node(Split(), name = 'split_masks') - split_masks.inputs.splits = [1, len(self.run_list) - 1] - split_masks.inputs.squeeze = True # Unfold one-element splits removing the list - subject_level.connect(select_files, 'masks', split_masks, 'inlist') + # Merge Node - Merge masks files for each subject + merge_masks = Node(Merge(), name = 'merge_masks') + merge_masks.inputs.dimension = 't' + subject_level.connect(select_files, 'masks', merge_masks, 'in_files') - # MultiImageMaths Node - Create a subject mask by + # MathsCommand Node - Create a subject mask by # computing the intersection of all run masks. - mask_intersection = Node(MultiImageMaths(), name = 'mask_intersection') - mask_intersection.inputs.op_string = '-mul %s ' * (len(self.run_list) - 1) - subject_level.connect(split_masks, 'out1', mask_intersection, 'in_file') - subject_level.connect(split_masks, 'out2', mask_intersection, 'operand_files') + mask_intersection = Node(MathsCommand(), name = 'mask_intersection') + mask_intersection.inputs.args = '-Tmin -thr 0.9' + subject_level.connect(merge_masks, 'merged_file', mask_intersection, 'in_file') # L2Model Node - Generate subject specific second level model generate_model = Node(L2Model(), name = 'generate_model') @@ -560,20 +558,17 @@ def get_group_level_analysis_sub_workflow(self, method): merge_varcopes.inputs.dimension = 't' group_level.connect(get_varcopes, ('out_list', clean_list), merge_varcopes, 'in_files') - # Split Node - Split mask list to serve them as inputs of the MultiImageMaths node. - split_masks = Node(Split(), name = 'split_masks') - split_masks.inputs.splits = [1, (len(self.subject_list) * len(self.run_list)) - 1] - split_masks.inputs.squeeze = True # Unfold one-element splits removing the list - group_level.connect(get_masks, ('out_list', clean_list), split_masks, 'inlist') + # Merge Node - Merge mask files + merge_masks = Node(Merge(), name = 'merge_masks') + merge_masks.inputs.dimension = 't' + group_level.connect(get_masks, ('out_list', clean_list), merge_masks, 'in_files') - # MultiImageMaths Node - Create a subject mask by + # MathsCommand Node - Create a subject mask by # computing the intersection of all run masks # (all voxels included in at least 80% of masks across all runs). - mask_intersection = Node(MultiImageMaths(), name = 'mask_intersection') - mask_intersection.inputs.op_string = '-add %s ' * (len(self.subject_list) - 1)\ - + f' -thr {len(self.subject_list) * len(self.run_list) * 0.8} -bin' - group_level.connect(split_masks, 'out1', mask_intersection, 'in_file') - group_level.connect(split_masks, 'out2', mask_intersection, 'operand_files') + mask_intersection = Node(MathsCommand(), name = 'mask_intersection') + mask_intersection.inputs.args = '-Tmin -thr 0.9' + group_level.connect(merge_masks, 'merged_file', mask_intersection, 'in_file') # MultipleRegressDesign Node - Specify model specify_model = Node(MultipleRegressDesign(), name = 'specify_model') From a525dc274cbf4364a0ac8129caef3f44bf89b10f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Fri, 28 Jun 2024 11:48:27 +0200 Subject: [PATCH 11/12] [TEST] updates for B23O --- narps_open/pipelines/team_B5I6.py | 9 +++------ tests/pipelines/test_team_B5I6.py | 27 ++++++++++++++------------- 2 files changed, 17 insertions(+), 19 deletions(-) diff --git a/narps_open/pipelines/team_B5I6.py b/narps_open/pipelines/team_B5I6.py index 2fe73a97..4d6cd317 100644 --- a/narps_open/pipelines/team_B5I6.py +++ b/narps_open/pipelines/team_B5I6.py @@ -714,10 +714,8 @@ def get_group_level_outputs(self): 'contrast_id': self.contrast_list, 'method': ['equalRange', 'equalIndifference'], 'file': [ - '_cluster0/zstat1_pval.nii.gz', # TODO : output for randomise - '_cluster0/zstat1_threshold.nii.gz', - '_cluster1/zstat2_pval.nii.gz', - '_cluster1/zstat2_threshold.nii.gz', + 'randomise_tfce_corrp_tstat1.nii.gz', + 'randomise_tfce_corrp_tstat2.nii.gz', 'tstat1.nii.gz', 'tstat2.nii.gz', 'zstat1.nii.gz', @@ -738,8 +736,7 @@ def get_group_level_outputs(self): parameters = { 'contrast_id': self.contrast_list, 'file': [ - '_cluster0/zstat1_pval.nii.gz', # TODO : output for randomise - '_cluster0/zstat1_threshold.nii.gz', + 'randomise_tfce_corrp_tstat1.nii.gz', 'tstat1.nii.gz', 'zstat1.nii.gz' ] diff --git a/tests/pipelines/test_team_B5I6.py b/tests/pipelines/test_team_B5I6.py index 956ee9ec..f80405da 100644 --- a/tests/pipelines/test_team_B5I6.py +++ b/tests/pipelines/test_team_B5I6.py @@ -52,11 +52,11 @@ def test_outputs(): # 1 - 1 subject outputs pipeline.subject_list = ['001'] - helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*1, 4*6*1, 8*4*2 + 4*4, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*1, 4*6*1, 6*6*2 + 3*2, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*4, 4*6*4, 8*4*2 + 4*4, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*4, 4*6*4, 6*6*2 + 3*2, 18]) @staticmethod @mark.unit_test @@ -64,9 +64,10 @@ def test_subject_information(): """ Test the get_subject_information method """ # Get test files - test_file = join(Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') - test_file_2 = join( - Configuration()['directories']['test_data'], 'pipelines', 'events_resp.tsv') + test_file = abspath(join( + Configuration()['directories']['test_data'], 'pipelines', 'events.tsv')) + test_file_2 = abspath(join( + Configuration()['directories']['test_data'], 'pipelines', 'events_resp.tsv')) # Prepare several scenarii info_missed = PipelineTeamB5I6.get_subject_information(test_file) @@ -120,11 +121,11 @@ def test_confounds_file_no_outliers(temporary_data_dir): """ Test the get_confounds_file method in the case with no outliers """ # Get input and reference output file - confounds_file = join( - Configuration()['directories']['test_data'], 'pipelines', 'confounds.tsv') - reference_file = join( + confounds_file = abspath(join( + Configuration()['directories']['test_data'], 'pipelines', 'confounds.tsv')) + reference_file = abspath(join( Configuration()['directories']['test_data'], - 'pipelines', 'team_B5I6', 'out_confounds_no_outliers.tsv') + 'pipelines', 'team_B5I6', 'out_confounds_no_outliers.tsv')) # Create new confounds file confounds_node = Node(Function( @@ -152,12 +153,12 @@ def test_confounds_file_outliers(temporary_data_dir): """ Test the get_confounds_file method in the case with outliers """ # Get input and reference output file - confounds_file = join( + confounds_file = abspath(join( Configuration()['directories']['test_data'], - 'pipelines', 'team_B5I6', 'confounds_with_outliers.tsv') - reference_file = join( + 'pipelines', 'team_B5I6', 'confounds_with_outliers.tsv')) + reference_file = abspath(join( Configuration()['directories']['test_data'], - 'pipelines', 'team_B5I6', 'out_confounds_outliers.tsv') + 'pipelines', 'team_B5I6', 'out_confounds_outliers.tsv')) # Create new confounds file confounds_node = Node(Function( From 951409c765827626fd0e1a5ce4dded8ab9829b2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= Date: Fri, 28 Jun 2024 12:04:40 +0200 Subject: [PATCH 12/12] Test update [skip ci] --- tests/pipelines/test_team_B5I6.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/pipelines/test_team_B5I6.py b/tests/pipelines/test_team_B5I6.py index f80405da..62ed6aec 100644 --- a/tests/pipelines/test_team_B5I6.py +++ b/tests/pipelines/test_team_B5I6.py @@ -52,11 +52,11 @@ def test_outputs(): # 1 - 1 subject outputs pipeline.subject_list = ['001'] - helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*1, 4*6*1, 6*6*2 + 3*2, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*1, 4*6*1, 6*6*2 + 3*6, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*4, 4*6*4, 6*6*2 + 3*2, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 4*6*4*4, 4*6*4, 6*6*2 + 3*6, 18]) @staticmethod @mark.unit_test