From ad8dde7e0dee1aadf9110f42c40635bc2acd7955 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= <117362283+bclenet@users.noreply.github.com> Date: Wed, 13 Mar 2024 16:59:14 +0100 Subject: [PATCH 1/8] UK24 reproduction (#179) * [UK24] writing parts of the preprocessing workflow * UK24 testable preprocessing * Preprocessing outputs * Subject_level_analysis * [TEST] adding tests for UK24 * [TEST] testing methods of UK24 pipelines * [TEST] compare_2d_arrays helper & temporary directory fixture * [LINT] linting pipeline tests * [TEST] temporary directory fixture * [DOC] using /work as volume mapping point inside nipype container + typos * Removing debug prints * [LINT] propagate use of fixture in test files [skip ci] --- narps_open/pipelines/__init__.py | 2 +- narps_open/pipelines/team_UK24.py | 886 ++++++++++++++++++ tests/pipelines/test_team_UK24.py | 207 ++++ .../team_UK24/average_values_csf.txt | 4 + .../pipelines/team_UK24/average_values_wm.txt | 4 + .../pipelines/team_UK24/confounds.tsv | 4 + .../team_UK24/framewise_displacement.txt | 4 + .../pipelines/team_UK24/func_resampled-32.nii | Bin 0 -> 2752 bytes .../pipelines/team_UK24/mask_resampled-32.nii | Bin 0 -> 1552 bytes .../team_UK24/realignment_parameters.txt | 4 + .../team_UK24/reference_average_values.txt | 4 + 11 files changed, 1118 insertions(+), 1 deletion(-) create mode 100755 narps_open/pipelines/team_UK24.py create mode 100644 tests/pipelines/test_team_UK24.py create mode 100644 tests/test_data/pipelines/team_UK24/average_values_csf.txt create mode 100644 tests/test_data/pipelines/team_UK24/average_values_wm.txt create mode 100644 tests/test_data/pipelines/team_UK24/confounds.tsv create mode 100644 tests/test_data/pipelines/team_UK24/framewise_displacement.txt create mode 100644 tests/test_data/pipelines/team_UK24/func_resampled-32.nii create mode 100644 tests/test_data/pipelines/team_UK24/mask_resampled-32.nii create mode 100644 tests/test_data/pipelines/team_UK24/realignment_parameters.txt create mode 100644 tests/test_data/pipelines/team_UK24/reference_average_values.txt diff --git a/narps_open/pipelines/__init__.py b/narps_open/pipelines/__init__.py index 2d67eeee..a9b62d9b 100644 --- a/narps_open/pipelines/__init__.py +++ b/narps_open/pipelines/__init__.py @@ -71,7 +71,7 @@ 'T54A': 'PipelineTeamT54A', 'U26C': 'PipelineTeamU26C', 'UI76': None, - 'UK24': None, + 'UK24': 'PipelineTeamUK24', 'V55J': None, 'VG39': None, 'X19V': 'PipelineTeamX19V', diff --git a/narps_open/pipelines/team_UK24.py b/narps_open/pipelines/team_UK24.py new file mode 100755 index 00000000..f4e8192f --- /dev/null +++ b/narps_open/pipelines/team_UK24.py @@ -0,0 +1,886 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Write the work of NARPS' team UK24 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.utility.base import Select, Merge, Split +from nipype.interfaces.io import SelectFiles, DataSink +from nipype.interfaces.spm import ( + Coregister, Segment, Reslice, Realign, + Smooth, Level1Design, OneSampleTTestDesign, TwoSampleTTestDesign, + EstimateModel, EstimateContrast, Threshold + ) +from nipype.algorithms.confounds import FramewiseDisplacement +from nipype.algorithms.modelgen import SpecifySPMModel +from nipype.algorithms.misc import Gunzip, SimpleThreshold + +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 ( + remove_file, list_intersection, elements_in_string, clean_list + ) + +class PipelineTeamUK24(Pipeline): + """ A class that defines the pipeline of team UK24. """ + + def __init__(self): + super().__init__() + self.fwhm = 4.0 + self.team_id = 'UK24' + self.contrast_list = ['0001', '0002'] + condition_names = [ + 'gain', 'gainxgain_value^1', 'gainxgain_rt^1', + 'loss', 'lossxloss_value^1', 'lossxloss_rt^1', + 'no_gain_no_loss', 'no_gain_no_lossxno_gain_no_loss_rt^1' + ] + self.subject_level_contrasts = [ + ['gain_param', 'T', condition_names, [0, 1, 1, 0, 0, 0, 0, 0]], + ['loss_param', 'T', condition_names, [0, 0, 0, 0, 1, 1, 0, 0]] + ] + + def get_average_values( + in_file: str, + mask: str, + subject_id: str, + run_id: str, + out_file_suffix: str + ) -> str: + """ + Compute the average value of each frame of a 4D image, while applying a mask. + Write the values in a text file. + + Parameters: + - in_file: str, path to the input file (4D) + - mask: str, path to the mask file (3D) Zeroes in mask are considered as 0, + other values are 1 + - subject_id: str, related subject id + - run_id: str, related run id + - out_file_suffix: str, string to append at the end of the file name + + Returns: + - preprocessing : nipype.WorkFlow + """ + from os import makedirs + from os.path import join, abspath + + from nilearn.image import load_img, math_img, index_img, get_data + + # Load input images + mask_img = load_img(mask) + in_file_img = load_img(in_file) + + # Loop through time points + average_values = [] + for frame_index in range(in_file_img.shape[3]): + average_values.append( + get_data(math_img( + 'np.mean(image * (mask > 0.0))', + image = index_img(in_file_img, frame_index), + mask = mask_img + )) + ) + + # Write confounds to a file + out_file_name = abspath(f'sub-{subject_id}_run-{run_id}_' + out_file_suffix) + + # Write output file + with open(out_file_name, 'w', encoding = 'utf-8') as writer: + for value in average_values: + writer.write(f'{value}\n') + + return out_file_name + + def get_preprocessing(self): + """ + Create the preprocessing workflow. + + Returns: + - preprocessing : nipype.WorkFlow + """ + # Initialize preprocessing workflow to connect nodes along the way + preprocessing = Workflow( + base_dir = self.directories.working_dir, + name = 'preprocessing') + + # IDENTITY INTERFACE - To iterate on subjects + information_source = Node(IdentityInterface( + fields = ['subject_id']), + name = 'information_source') + information_source.iterables = [('subject_id', self.subject_list)] + + # SELECT FILES - to select necessary files + templates = { + 'anat' : join('sub-{subject_id}', 'anat', 'sub-{subject_id}_T1w.nii.gz'), + 'func' : join('sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-*_bold.nii.gz'), + 'sbref' : join('sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-*_sbref.nii.gz'), + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory = self.directories.dataset_dir + preprocessing.connect(information_source, 'subject_id', select_files, 'subject_id') + + # GUNZIP - gunzip files because SPM do not use .nii.gz files + gunzip_anat = Node(Gunzip(), name = 'gunzip_anat') + gunzip_func = MapNode(Gunzip(), name = 'gunzip_func', iterfield = ['in_file']) + gunzip_sbref = MapNode(Gunzip(), name = 'gunzip_sbref', iterfield = ['in_file']) + preprocessing.connect(select_files, 'anat', gunzip_anat, 'in_file') + preprocessing.connect(select_files, 'func', gunzip_func, 'in_file') + preprocessing.connect(select_files, 'sbref', gunzip_sbref, 'in_file') + + # SELECT - Select the first sbref file from the selected (and gunzipped) files. + select_first_sbref = Node(Select(), name = 'select_first_sbref') + select_first_sbref.inputs.index = [0] + preprocessing.connect(gunzip_sbref, 'out_file', select_first_sbref, 'inlist') + + # COREGISTER - Coregistration of the structural T1w image to the reference functional image + # (defined as the single volume 'sbref' image acquired before the first functional run). + coregistration = Node(Coregister(), name = 'coregistration') + coregistration.inputs.cost_function = 'nmi' + preprocessing.connect(gunzip_anat, 'out_file', coregistration, 'source') + preprocessing.connect(select_first_sbref, 'out', coregistration, 'target') + + # SEGMENT - Segmentation of the coregistered T1w image into grey matter, white matter and + # cerebrospinal fluid tissue maps. + segmentation = Node(Segment(), name = 'segmentation') + segmentation.inputs.bias_fwhm = 60 #mm + #segmentation.inputs.bias_regularization = (0 or 1e-05 or 0.0001 or 0.001 or 0.01 or 0.1 or 1 or 10) "light regularization" + segmentation.inputs.csf_output_type = [False,False,True] # Output maps in native space only + segmentation.inputs.gm_output_type = [False,False,True] # Output maps in native space only + segmentation.inputs.wm_output_type = [False,False,True] # Output maps in native space only + preprocessing.connect(coregistration, 'coregistered_source', segmentation, 'data') + + # MERGE - Merge files for the reslicing node into one input. + merge_before_reslicing = Node(Merge(4), name = 'merge_before_reslicing') + preprocessing.connect(coregistration, 'coregistered_source', merge_before_reslicing, 'in1') + preprocessing.connect(segmentation, 'native_csf_image', merge_before_reslicing, 'in2') + preprocessing.connect(segmentation, 'native_wm_image', merge_before_reslicing, 'in3') + preprocessing.connect(segmentation, 'native_gm_image', merge_before_reslicing, 'in4') + + # RESLICE - Reslicing of coregistered T1w image and segmentations to the same voxel space + # of the reference (sbref) image. + reslicing = MapNode(Reslice(), name = 'reslicing', iterfield = 'in_file') + preprocessing.connect(merge_before_reslicing, 'out', reslicing, 'in_file') + preprocessing.connect(select_first_sbref, 'out', reslicing, 'space_defining') + + # REALIGN - Realigning all 4 sbref files images to the one acquired before the first run + # Realigns to the first image of the list by default. TODO : check + # Note : gunzip_sbref.out_file is a list of files + realign_sbref = Node(Realign(), name = 'realign_sbref') + realign_sbref.inputs.register_to_mean = False + preprocessing.connect(gunzip_sbref, 'out_file', realign_sbref, 'in_files') + + # REALIGN - Realigning all volumes in a particular run to the first image in that run. + # Realigns to the first image of the run by default. + # Note : gunzip_func.out_file is a list of files, but we want to realign run by run, + # therefore, passing one func file at a time. + realign_func = MapNode(Realign(), name = 'realign_func', iterfield = 'in_files') + realign_func.inputs.register_to_mean = False + preprocessing.connect(gunzip_func, 'out_file', realign_func, 'in_files') + + # SMOOTH - Spatial smoothing of fMRI data. + # Note : realign_func.realigned_files will be a list(list(files)) : + # we need a MapNode to process it. + smoothing = MapNode(Smooth(), name = 'smoothing', iterfield = 'in_files') + smoothing.inputs.fwhm = [self.fwhm] * 3 + preprocessing.connect(realign_func, 'realigned_files', smoothing, 'in_files') + + # SELECT - Select the segmentation probability maps. + select_maps = Node(Select(), name = 'select_maps') + select_maps.inputs.index = [1, 2] # csf and wm only + preprocessing.connect(reslicing, 'out_file', select_maps, 'inlist') + + # SIMPLE THRESHOLD - Apply a threshold to proability maps + threshold = Node(SimpleThreshold(), name = 'threshold') + threshold.inputs.threshold = 0.5 + preprocessing.connect(select_maps, 'out', threshold, 'volumes') + + # FRAMEWISE DISPLACEMENT - Compute framewise displacement from realignment parameters. + displacement = MapNode( + FramewiseDisplacement(), + name = 'displacement', + iterfield = 'in_file' + ) + displacement.inputs.parameter_source = 'SPM' + preprocessing.connect(realign_func, 'realignment_parameters', displacement, 'in_file') + + # SPLIT - Select the WM and CSF masks. + split_masks = Node(Split(), name = 'split_masks') + split_masks.inputs.splits = [1, 1] + split_masks.inputs.squeeze = True # Unfold one-element splits removing the list. + preprocessing.connect(threshold, 'thresholded_volumes', split_masks, 'inlist') + + # FUNCTION - Create text files with average values in WM at each timepoint + # to be later used as regressor. + average_values_wm = MapNode( + Function( + function = self.get_average_values, + input_names = ['in_file', 'mask', 'subject_id', 'run_id', 'out_file_suffix'], + output_names = ['out_file'] + ), + iterfield = ['run_id', 'in_file'], + name = 'average_values_wm') + average_values_wm.inputs.run_id = self.run_list + average_values_wm.inputs.out_file_suffix = join('reg_wm.txt') + preprocessing.connect(information_source, 'subject_id', average_values_wm, 'subject_id') + preprocessing.connect(smoothing, 'smoothed_files', average_values_wm, 'in_file') + preprocessing.connect(split_masks, 'out2', average_values_wm, 'mask') + + # FUNCTION - Create text files with average values in CSF at each timepoint + # to be later used as regressor. + average_values_csf = MapNode( + Function( + function = self.get_average_values, + input_names = ['in_file', 'mask', 'subject_id', 'run_id', 'out_file_suffix'], + output_names = ['out_file'] + ), + iterfield = ['run_id', 'in_file'], + name = 'average_values_csf') + average_values_csf.inputs.run_id = self.run_list + average_values_csf.inputs.out_file_suffix = join('reg_csf.txt') + preprocessing.connect(information_source, 'subject_id', average_values_csf, 'subject_id') + preprocessing.connect(smoothing, 'smoothed_files', average_values_csf, 'in_file') + preprocessing.connect(split_masks, 'out1', average_values_csf, 'mask') + + # DATA SINK - store the wanted results in the wanted repository + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + preprocessing.connect(smoothing, 'smoothed_files', data_sink, 'preprocessing.@smoothed') + preprocessing.connect(segmentation, 'transformation_mat', data_sink, 'preprocessing.@tf_to_mni') + preprocessing.connect(displacement, 'out_file', data_sink, 'preprocessing.@reg_scrubbing') + preprocessing.connect( + realign_func, 'realignment_parameters', + data_sink, 'preprocessing.@realignment_parameters') + preprocessing.connect(average_values_csf, 'out_file', data_sink, 'preprocessing.@reg_csf') + preprocessing.connect(average_values_wm, 'out_file', data_sink, 'preprocessing.@reg_wm') + + return preprocessing + + def get_preprocessing_outputs(self): + """ Return the names of the files the preprocessing is supposed to generate. """ + + output_dir = join(self.directories.output_dir, 'preprocessing', '_subject_id_{subject_id}') + + # WM and CSF average values + templates = [join(output_dir, f'_average_values_wm{index}', + 'sub-{subject_id}'+f'_run-{run_id}_reg_wm.txt')\ + for index, run_id in zip(range(len(self.run_list)), self.run_list)] + templates += [join(output_dir, f'_average_values_csf{index}', + 'sub-{subject_id}'+f'_run-{run_id}_reg_csf.txt')\ + for index, run_id in zip(range(len(self.run_list)), self.run_list)] + + # Framewise displacement values + templates += [join(output_dir, f'_displacement{index}', + 'fd_power_2012.txt')\ + for index in range(len(self.run_list))] + + # Realignement parameters + templates += [join(output_dir, f'_realign_func{index}', + 'rp_sub-{subject_id}'+f'_task-MGT_run-{run_id}_bold.txt')\ + for index, run_id in zip(range(len(self.run_list)), self.run_list)] + + # Smoothing outputs + templates += [join(output_dir, f'_smoothing{index}', + 'srsub-{subject_id}'+f'_task-MGT_run-{run_id}_bold.nii')\ + for index, run_id in zip(range(len(self.run_list)), self.run_list)] + + + # Transformation matrix to MNI + templates.append(join(output_dir, 'rsub-{subject_id}_T1w_seg_sn.mat')) + + # 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_run_level_analysis(self): + """ No run level analysis has been done by team UK24 """ + return None + + def get_subject_information(event_file): + """ + Create Bunchs for specifySPMModel, from data extracted from an event_file. + + Parameters : + - event_files: str, event file (one per run) for the subject + + Returns : + - subject_information: Bunch, relevant event information for subject level analysis. + """ + from nipype.interfaces.base import Bunch + from nipype.algorithms.modelgen import orth + + # Init empty lists inside directiries + onsets = [] + durations = [] + onsets_no_gain_no_loss = [] + durations_no_gain_no_loss = [] + gain_value = [] + gain_rt = [] + loss_value = [] + loss_rt = [] + no_gain_no_loss_rt = [] + + with open(event_file, 'rt') as file: + next(file) # skip the header + + for line in file: + info = line.strip().split() + + if 'accept' in info[5]: + onsets.append(float(info[0])) + durations.append(float(info[1])) + gain_value.append(float(info[2])) + loss_value.append(float(info[3])) + gain_rt.append(float(info[4])) + loss_rt.append(float(info[4])) + else: + onsets_no_gain_no_loss.append(float(info[0])) + durations_no_gain_no_loss.append(float(info[1])) + no_gain_no_loss_rt.append(float(info[4])) + + # NOTE : SPM automatically mean-centers regressors + #gain_value = list(gain_value - mean(gain_value)) + #loss_value = list(loss_value - mean(loss_value)) + #gain_rt = list(gain_rt - mean(gain_rt)) + #loss_rt = list(loss_rt - mean(loss_rt)) + #no_gain_no_loss_rt = list(no_gain_no_loss_rt - mean(no_gain_no_loss_rt)) + + # Othogonalize + #gain_rt = orth(gain_value, gain_rt) + #loss_rt = orth(loss_value, loss_rt) + + # Fill Bunch + return Bunch( + conditions = ['gain', 'loss', 'no_gain_no_loss'], + onsets = [onsets, onsets, onsets_no_gain_no_loss], + durations = [durations, durations, durations_no_gain_no_loss], + amplitudes = None, + tmod = None, + pmod = [ + Bunch( + name = ['gain_value', 'gain_rt'], + poly = [1, 1], + param = [gain_value, gain_rt] + ), + Bunch( + name = ['loss_value', 'loss_rt'], + poly = [1, 1], + param = [loss_value, loss_rt] + ), + Bunch( + name = ['no_gain_no_loss_rt'], + poly = [1], + param = [no_gain_no_loss_rt] + ) + ], + regressor_names = None, + regressors = None + ) + + def get_confounds_file( + framewise_displacement_file: str, + wm_average_file: str, + csf_average_file: str, + realignement_parameters: str, + subject_id: str, + run_id: str) -> str: + """ + Create a tsv file with only desired confounds per subject per run. + + Parameters : + - framewise_displacement_file: str, path to the framewise displacement file + - wm_average_file: str, path to the WM average values file + - csf_average_file: str, path to the CSF average values file + - realignement_parameters : path to the realignment parameters file + - subject_id : related subject id + - run_id : related run id + + Return : + - confounds_file : path to new file containing only desired confounds + """ + from os.path import abspath + + from pandas import DataFrame, read_csv + from numpy import array, transpose, concatenate, insert + + # Get the dataframe containing the 6 head motion parameter regressors + realign_data_frame = array(read_csv(realignement_parameters, sep = r'\s+', header = None)) + + # Get the dataframes containing the 2 tissue signal regressors + regressor_csf_data_frame = array(read_csv(csf_average_file, sep = '\t', header = None)) + regressor_wm_data_frame = array(read_csv(wm_average_file, sep = '\t', header = None)) + + # Get the dataframe containing framewise displacement + # and transform it as a scrubbing regressor. + scrubbing_data_frame = insert( + array( + (read_csv(framewise_displacement_file, sep = '\t', header = 0) > 0.5).astype(int) + ), + 0, 0, axis = 0) # Add a value of 0 at the beginning (first frame) + + # Extract all parameters + retained_parameters = DataFrame( + concatenate( + (realign_data_frame, + regressor_csf_data_frame, regressor_wm_data_frame, scrubbing_data_frame), + axis = 1) + ) + + # 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_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 + WARNING: the name attribute of the workflow is 'subject_level_analysis' + """ + # Create subject level analysis workflow + subject_level = Workflow( + base_dir = self.directories.working_dir, + name = 'subject_level_analysis') + + # IDENTITY INTERFACE - To iterate on subjects + information_source = Node(IdentityInterface( + fields = ['subject_id']), + name = 'information_source') + information_source.iterables = [('subject_id', self.subject_list)] + + # SELECT FILES - to select necessary files + templates = { + 'framewise_displacement_files' : join(self.directories.output_dir, 'preprocessing', + '_subject_id_{subject_id}', '_displacement*', + 'fd_power_2012.txt'), + 'wm_average_files' : join(self.directories.output_dir, 'preprocessing', + '_subject_id_{subject_id}', '_average_values_wm*', + 'sub-{subject_id}_run-*_reg_wm.txt'), + 'csf_average_files' : join(self.directories.output_dir, 'preprocessing', + '_subject_id_{subject_id}', '_average_values_csf*', + 'sub-{subject_id}_run-*_reg_csf.txt'), + 'realignement_parameters' : join(self.directories.output_dir, 'preprocessing', + '_subject_id_{subject_id}', '_realign_func*', + 'rp_sub-{subject_id}_task-MGT_run-*_bold.txt'), + 'func' : join(self.directories.output_dir, 'preprocessing', '_subject_id_{subject_id}', + '_smoothing*', 'srsub-{subject_id}_task-MGT_run-*_bold.nii'), + 'event' : 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 + subject_level.connect(information_source, 'subject_id', select_files, 'subject_id') + + # FUNCTION get_subject_information - generate files with event data + subject_information = MapNode(Function( + function = self.get_subject_information, + input_names = ['event_file'], + output_names = ['subject_info']), + iterfield = 'event_file', + name = 'subject_information') + subject_level.connect(select_files, 'event', subject_information, 'event_file') + + # FUNCTION node get_confounds_file - generate files with confounds data + confounds = MapNode( + Function( + function = self.get_confounds_file, + input_names = ['framewise_displacement_file', 'wm_average_file', + 'csf_average_file', 'realignement_parameters', 'subject_id', 'run_id'], + output_names = ['confounds_file'] + ), + name = 'confounds', + iterfield = ['framewise_displacement_file', 'wm_average_file', + 'csf_average_file', 'realignement_parameters', 'run_id']) + confounds.inputs.run_id = self.run_list + subject_level.connect(information_source, 'subject_id', confounds, 'subject_id') + subject_level.connect( + select_files, 'framewise_displacement_files', confounds, 'framewise_displacement_file') + subject_level.connect(select_files, 'wm_average_files', confounds, 'wm_average_file') + subject_level.connect(select_files, 'csf_average_files', confounds, 'csf_average_file') + subject_level.connect( + select_files, 'realignement_parameters', confounds, 'realignement_parameters') + + # SPECIFY MODEL - generates SPM-specific Model + specify_model = Node(SpecifySPMModel(), name = 'specify_model') + #specify_model.inputs.concatenate_runs = True # TODO : actually concatenate runs ???? + 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 + subject_level.connect(select_files, 'func', specify_model, 'functional_runs') + subject_level.connect(confounds, 'confounds_file', specify_model, 'realignment_parameters') + subject_level.connect(subject_information, 'subject_info', specify_model, 'subject_info') + + # LEVEL 1 DESIGN - Generates an SPM design matrix + model_design = Node(Level1Design(), name = 'model_design') + model_design.inputs.bases = {'hrf': {'derivs': [0, 0]}} + model_design.inputs.timing_units = 'secs' + model_design.inputs.interscan_interval = TaskInformation()['RepetitionTime'] + subject_level.connect(specify_model, 'session_info', model_design, 'session_info') + + # ESTIMATE MODEL - estimate the parameters of the model + model_estimate = Node(EstimateModel(), name = 'model_estimate') + model_estimate.inputs.estimation_method = {'Classical': 1} + subject_level.connect(model_design, 'spm_mat_file', model_estimate, 'spm_mat_file') + + # ESTIMATE CONTRAST - estimates contrasts + contrast_estimate = Node(EstimateContrast(), name = 'contrast_estimate') + contrast_estimate.inputs.contrasts = self.subject_level_contrasts + subject_level.connect([ + (model_estimate, contrast_estimate, [ + ('spm_mat_file', 'spm_mat_file'), + ('beta_images', 'beta_images'), + ('residual_image', 'residual_image') + ]) + ]) + + # DATA SINK - 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, data_sink, [ + ('con_images', 'subject_level_analysis.@con_images'), + ('spmT_images', 'subject_level_analysis.@spmT_images'), + ('spm_mat_file', 'subject_level_analysis.@spm_mat_file') + ]) + ]) + + return subject_level + + def get_subject_level_outputs(self): + """ Return the names of the files the subject level analysis is supposed to generate. """ + + # Contrat maps + templates = [join( + self.directories.output_dir, + 'subject_level_analysis', '_subject_id_{subject_id}', f'con_{contrast_id}.nii')\ + for contrast_id in self.contrast_list] + + # SPM.mat file + templates += [join( + self.directories.output_dir, + 'subject_level_analysis', '_subject_id_{subject_id}', 'SPM.mat')] + + # spmT maps + templates += [join( + self.directories.output_dir, + 'subject_level_analysis', '_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 + """ + + 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_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 + information_source = Node( + IdentityInterface( + fields=['contrast_id']), + name='information_source') + information_source.iterables = [('contrast_id', self.contrast_list)] + + # SelectFiles + templates = { + # Contrasts for all participants + 'contrasts' : join(self.directories.output_dir, + 'subject_level_analysis', '_subject_id_*', 'con_{contrast_id}.nii') + } + + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory = self.directories.results_dir + select_files.inputs.force_lists = True + + # Datasink - save important files + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + + # 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 + + # 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' + ) + + # Estimate model + estimate_model = Node(EstimateModel(), name = 'estimate_model') + estimate_model.inputs.estimation_method = {'Classical':1} + + # Estimate contrasts + estimate_contrast = Node(EstimateContrast(), name = 'estimate_contrast') + estimate_contrast.inputs.group_contrast = True + + ## Create thresholded maps + threshold = MapNode(Threshold(), + name = 'threshold', + iterfield = ['stat_image', 'contrast_index']) + threshold.inputs.use_fwe_correction = False + threshold.inputs.height_threshold = 0.001 + threshold.inputs.extent_fdr_p_threshold = 0.05 + threshold.inputs.use_topo_fdr = False + threshold.inputs.force_activation = True + threshold.synchronize = True + + group_level_analysis = Workflow( + base_dir = self.directories.working_dir, + name = f'group_level_analysis_{method}_nsub_{nb_subjects}') + group_level_analysis.connect([ + (information_source, select_files, [('contrast_id', 'contrast_id')]), + (select_files, get_contrasts, [('contrasts', 'input_str')]), + (estimate_model, estimate_contrast, [ + ('spm_mat_file', 'spm_mat_file'), + ('residual_image', 'residual_image'), + ('beta_images', 'beta_images')]), + (estimate_contrast, threshold, [ + ('spm_mat_file', 'spm_mat_file'), + ('spmT_images', 'stat_image')]), + (estimate_model, data_sink, [ + ('mask_image', f'group_level_analysis_{method}_nsub_{nb_subjects}.@mask')]), + (estimate_contrast, data_sink, [ + ('spm_mat_file', f'group_level_analysis_{method}_nsub_{nb_subjects}.@spm_mat'), + ('spmT_images', f'group_level_analysis_{method}_nsub_{nb_subjects}.@T'), + ('con_images', f'group_level_analysis_{method}_nsub_{nb_subjects}.@con')]), + (threshold, data_sink, [ + ('thresholded_map', f'group_level_analysis_{method}_nsub_{nb_subjects}.@thresh')])]) + + if method in ('equalRange', 'equalIndifference'): + estimate_contrast.inputs.contrasts = [ + ('Group', 'T', ['mean'], [1]), + ('Group', 'T', ['mean'], [-1]) + ] + + threshold.inputs.contrast_index = [1, 2] + + # Specify design matrix + one_sample_t_test_design = Node(OneSampleTTestDesign(), + name = 'one_sample_t_test_design') + group_level_analysis.connect([ + (get_contrasts, one_sample_t_test_design, [ + (('out_list', clean_list), 'in_files') + ]), + (one_sample_t_test_design, estimate_model, [('spm_mat_file', 'spm_mat_file')]) + ]) + + if method == 'equalRange': + group_level_analysis.connect([ + (get_equal_range_subjects, get_contrasts, [ + (('out_list', complete_subject_ids), 'elements') + ]) + ]) + + elif method == 'equalIndifference': + group_level_analysis.connect([ + (get_equal_indifference_subjects, get_contrasts, [ + (('out_list', complete_subject_ids), 'elements') + ]) + ]) + + elif method == 'groupComp': + estimate_contrast.inputs.contrasts = [ + ('Eq range vs Eq indiff in loss', 'T', ['Group_{1}', 'Group_{2}'], [1, -1]) + ] + + threshold.inputs.contrast_index = [1] + + # 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_2 = MapNode(Function( + function = elements_in_string, + input_names = ['input_str', 'elements'], + output_names = ['out_list'] + ), + name = 'get_contrasts_2', iterfield = 'input_str' + ) + + # Specify design matrix + two_sample_t_test_design = Node(TwoSampleTTestDesign(), + name = 'two_sample_t_test_design') + two_sample_t_test_design.inputs.unequal_variance = True + + group_level_analysis.connect([ + (select_files, get_contrasts_2, [('contrasts', 'input_str')]), + (get_equal_range_subjects, get_contrasts, [ + (('out_list', complete_subject_ids), 'elements') + ]), + (get_equal_indifference_subjects, get_contrasts_2, [ + (('out_list', complete_subject_ids), 'elements') + ]), + (get_contrasts, two_sample_t_test_design, [ + (('out_list', clean_list), 'group1_files') + ]), + (get_contrasts_2, two_sample_t_test_design, [ + (('out_list', clean_list), 'group2_files') + ]), + (two_sample_t_test_design, estimate_model, [('spm_mat_file', 'spm_mat_file')]) + ]) + + 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. + Note that hypotheses 5 to 8 correspond to the maps given by the team in their results ; + but they are not fully consistent with the hypotheses definitions as expected by NARPS. + """ + nb_sub = len(self.subject_list) + files = [ + 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'), + 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'), + 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'), + 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'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0003', '_threshold1', 'spmT_0002_thr.nii'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0003', 'spmT_0002.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0003', '_threshold1', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0003', 'spmT_0001.nii'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0003', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0003', 'spmT_0001.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0003', '_threshold0', 'spmT_0002_thr.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0003', 'spmT_0002.nii'), + join(f'group_level_analysis_groupComp_nsub_{nb_sub}', + '_contrast_id_0003', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_groupComp_nsub_{nb_sub}', + '_contrast_id_0003', 'spmT_0001.nii') + ] + return [join(self.directories.output_dir, f) for f in files] diff --git a/tests/pipelines/test_team_UK24.py b/tests/pipelines/test_team_UK24.py new file mode 100644 index 00000000..7d420c02 --- /dev/null +++ b/tests/pipelines/test_team_UK24.py @@ -0,0 +1,207 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Tests of the 'narps_open.pipelines.team_UK24' module. + +Launch this test with PyTest + +Usage: +====== + pytest -q test_team_UK24.py + pytest -q test_team_UK24.py -k +""" + +from os.path import join, exists, abspath +from filecmp import cmp + +from pytest import helpers, mark +from nipype import Workflow, Node +from nipype.interfaces.utility import Function +from nipype.interfaces.base import Bunch + +from narps_open.utils.configuration import Configuration +from narps_open.pipelines.team_UK24 import PipelineTeamUK24 + +class TestPipelinesTeamUK24: + """ A class that contains all the unit tests for the PipelineTeamUK24 class.""" + + @staticmethod + @mark.unit_test + def test_create(): + """ Test the creation of a PipelineTeamUK24 object """ + + pipeline = PipelineTeamUK24() + + # 1 - check the parameters + assert pipeline.fwhm == 4.0 + assert pipeline.team_id == 'UK24' + + # 2 - check workflows + assert isinstance(pipeline.get_preprocessing(), Workflow) + 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 PipelineTeamUK24 object """ + pipeline = PipelineTeamUK24() + # 1 - 1 subject outputs + pipeline.subject_list = ['001'] + helpers.test_pipeline_outputs(pipeline, [21, 0, 5, 42, 18]) + + # 2 - 4 subjects outputs + pipeline.subject_list = ['001', '002', '003', '004'] + helpers.test_pipeline_outputs(pipeline, [84, 0, 20, 42, 18]) + + @staticmethod + @mark.unit_test + def test_average_values(temporary_data_dir): + """ Test the get_average_values method + For this test, we created the two following files by downsampling output files + from preprocessing nodes : + - func_resampled-32.nii (smoothed functional files) + - mask_resampled-32.nii (white matter mask) + Voxel dimension was multiplied by 32, number of slices was reduced to 4. + """ + # Test files + test_in_file = abspath(join(Configuration()['directories']['test_data'], + 'pipelines', 'team_UK24', 'func_resampled-32.nii')) + test_mask_file = abspath(join(Configuration()['directories']['test_data'], + 'pipelines', 'team_UK24', 'mask_resampled-32.nii')) + + # Reference file + ref_file = abspath(join(Configuration()['directories']['test_data'], + 'pipelines', 'team_UK24', 'reference_average_values.txt')) + + # Create average values file + average_node = Node(Function( + input_names = ['in_file', 'mask', 'subject_id', 'run_id', 'out_file_suffix'], + output_names = ['out_file'], + function = PipelineTeamUK24.get_average_values), + name = 'average_node') + average_node.base_dir = temporary_data_dir + average_node.inputs.in_file = test_in_file + average_node.inputs.mask = test_mask_file + average_node.inputs.subject_id = 'sid' + average_node.inputs.run_id = 'rid' + average_node.inputs.out_file_suffix = 'suffix.txt' + average_node.run() + + # Check file was created + created_average_values_file = abspath(join( + temporary_data_dir, average_node.name, 'sub-sid_run-rid_suffix.txt')) + assert exists(created_average_values_file) + + # Check contents + assert cmp(ref_file, created_average_values_file) + + @staticmethod + @mark.unit_test + def test_confounds_file(temporary_data_dir): + """ Test the get_confounds_file method """ + + # Test files + test_average_values_csf = abspath(join(Configuration()['directories']['test_data'], + 'pipelines', 'team_UK24', 'average_values_csf.txt')) + test_average_values_wm = abspath(join(Configuration()['directories']['test_data'], + 'pipelines', 'team_UK24', 'average_values_wm.txt')) + test_framewise_displacement = abspath(join(Configuration()['directories']['test_data'], + 'pipelines', 'team_UK24', 'framewise_displacement.txt')) + test_realignment_parameters = abspath(join(Configuration()['directories']['test_data'], + 'pipelines', 'team_UK24', 'realignment_parameters.txt')) + + # Reference file + ref_file = abspath(join(Configuration()['directories']['test_data'], + 'pipelines', 'team_UK24', 'confounds.tsv')) + + # Create average values file + confounds_node = Node(Function( + input_names = [ + 'framewise_displacement_file', 'wm_average_file', 'csf_average_file', + 'realignement_parameters', 'subject_id', 'run_id' + ], + output_names = ['out_file'], + function = PipelineTeamUK24.get_confounds_file), + name = 'confounds_node') + confounds_node.base_dir = temporary_data_dir + confounds_node.inputs.framewise_displacement_file = test_framewise_displacement + confounds_node.inputs.wm_average_file = test_average_values_wm + confounds_node.inputs.csf_average_file = test_average_values_csf + confounds_node.inputs.realignement_parameters = test_realignment_parameters + confounds_node.inputs.subject_id = 'sid' + confounds_node.inputs.run_id = 'rid' + confounds_node.run() + + # Check 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) + + # Check contents + assert cmp(ref_file, created_confounds_file) + + @staticmethod + @mark.unit_test + def test_subject_information(): + """ Test the get_subject_information method """ + + # Test with 'gain' + test_event_file = join( + Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') + information = PipelineTeamUK24.get_subject_information(test_event_file) + + assert isinstance(information, Bunch) + assert information.conditions == ['gain', 'loss', 'no_gain_no_loss'] + + helpers.compare_float_2d_arrays([ + [4.0, 4.0], + [4.0, 4.0], + [4.0, 4.0, 4.0] + ], + information.durations) + helpers.compare_float_2d_arrays([ + [4.071, 11.834], + [4.071, 11.834], + [19.535, 27.535, 36.435] + ], + information.onsets) + + paramateric_modulation = information.pmod[0] + assert isinstance(paramateric_modulation, Bunch) + assert paramateric_modulation.name == ['gain_value', 'gain_rt'] + assert paramateric_modulation.poly == [1, 1] + helpers.compare_float_2d_arrays([ + [14.0, 34.0], + [2.388, 2.289] + ], + paramateric_modulation.param) + + paramateric_modulation = information.pmod[1] + assert isinstance(paramateric_modulation, Bunch) + assert paramateric_modulation.name == ['loss_value', 'loss_rt'] + assert paramateric_modulation.poly == [1, 1] + helpers.compare_float_2d_arrays([ + [6.0, 14.0], + [2.388, 2.289] + ], + paramateric_modulation.param) + + paramateric_modulation = information.pmod[2] + assert isinstance(paramateric_modulation, Bunch) + assert paramateric_modulation.name == ['no_gain_no_loss_rt'] + assert paramateric_modulation.poly == [1] + helpers.compare_float_2d_arrays([ + [0.0, 2.08, 2.288] + ], + paramateric_modulation.param) + + @staticmethod + @mark.pipeline_test + def test_execution(): + """ Test the execution of a PipelineTeamUK24 and compare results """ + helpers.test_pipeline_evaluation('UK24') diff --git a/tests/test_data/pipelines/team_UK24/average_values_csf.txt b/tests/test_data/pipelines/team_UK24/average_values_csf.txt new file mode 100644 index 00000000..6ee5b66b --- /dev/null +++ b/tests/test_data/pipelines/team_UK24/average_values_csf.txt @@ -0,0 +1,4 @@ +309.79840087890625 +306.74420166015625 +305.7932434082031 +305.5064392089844 diff --git a/tests/test_data/pipelines/team_UK24/average_values_wm.txt b/tests/test_data/pipelines/team_UK24/average_values_wm.txt new file mode 100644 index 00000000..c6fff4a5 --- /dev/null +++ b/tests/test_data/pipelines/team_UK24/average_values_wm.txt @@ -0,0 +1,4 @@ +618.2028198242188 +617.9356689453125 +618.0912475585938 +617.5615234375 diff --git a/tests/test_data/pipelines/team_UK24/confounds.tsv b/tests/test_data/pipelines/team_UK24/confounds.tsv new file mode 100644 index 00000000..4f4568df --- /dev/null +++ b/tests/test_data/pipelines/team_UK24/confounds.tsv @@ -0,0 +1,4 @@ +0.0 0.0 0.0 0.0 0.0 0.0 309.79840087890625 618.2028198242188 0.0 +0.002151997 -0.033891228 0.039675607 0.0013263004 1.1828492e-05 0.00034149533 306.74420166015625 617.9356689453125 0.0 +0.0004410012 -0.050145525 0.12030211 0.0015248415 0.00011492542 0.00069859071 305.7932434082031 618.0912475585938 0.0 +0.0025807056 0.0015129009 0.16351441 0.0012306972 0.00024110722 0.00056041322 305.5064392089844 617.5615234375 1.0 diff --git a/tests/test_data/pipelines/team_UK24/framewise_displacement.txt b/tests/test_data/pipelines/team_UK24/framewise_displacement.txt new file mode 100644 index 00000000..0239c25a --- /dev/null +++ b/tests/test_data/pipelines/team_UK24/framewise_displacement.txt @@ -0,0 +1,4 @@ +FramewiseDisplacement +1.597000431000000220e-01 +1.315284661999999993e-01 +8.017163291999997166e-01 diff --git a/tests/test_data/pipelines/team_UK24/func_resampled-32.nii b/tests/test_data/pipelines/team_UK24/func_resampled-32.nii new file mode 100644 index 0000000000000000000000000000000000000000..7958202f56af5e2e2e706b989daf3caadeafddb7 GIT binary patch literal 2752 zcma!HWFQK#Ft9SP0Wk{$BN(D;5@1jO3pCg>Ff=%U@T{3LK{SYugf}=t_^7%Ej>m*z zwbkEChtgfY9%}Zwd%i8b5eVym|WE@re$%C+)_ zlS}jUNVpg!E=&{RarK?Z>ef_%NH1}R_1vzze0QCCcr)A%EeTFnGx18dW7-{XF-q(m z)k~ic@cr-S`nERP?R_Ppykf3?>hkW|N|)Qm5ardwtrD(RW87T!Y)Ga`{05b9x?XwIfI79fTx(ANOgkrVT z-%E$mUB4b`_PTrs ML!)3c1ZW=u0HIkSJpcdz literal 0 HcmV?d00001 diff --git a/tests/test_data/pipelines/team_UK24/realignment_parameters.txt b/tests/test_data/pipelines/team_UK24/realignment_parameters.txt new file mode 100644 index 00000000..bafc627e --- /dev/null +++ b/tests/test_data/pipelines/team_UK24/realignment_parameters.txt @@ -0,0 +1,4 @@ + 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 + 2.1519970e-03 -3.3891228e-02 3.9675607e-02 1.3263004e-03 1.1828492e-05 3.4149533e-04 + 4.4100120e-04 -5.0145525e-02 1.2030211e-01 1.5248415e-03 1.1492542e-04 6.9859071e-04 + 2.5807056e-03 1.5129009e-03 1.6351441e-01 1.2306972e-03 2.4110722e-04 5.6041322e-04 diff --git a/tests/test_data/pipelines/team_UK24/reference_average_values.txt b/tests/test_data/pipelines/team_UK24/reference_average_values.txt new file mode 100644 index 00000000..f35a6487 --- /dev/null +++ b/tests/test_data/pipelines/team_UK24/reference_average_values.txt @@ -0,0 +1,4 @@ +52.34728240966797 +52.014034271240234 +52.75069808959961 +52.45905303955078 From 273a8d7c3a24ac67831de6d7fc42f228c9b6bf95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= <117362283+bclenet@users.noreply.github.com> Date: Wed, 13 Mar 2024 17:20:02 +0100 Subject: [PATCH 2/8] Correcting issue with wrong test parametrization (#184) --- tests/pipelines/test_team_98BT.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/tests/pipelines/test_team_98BT.py b/tests/pipelines/test_team_98BT.py index 6e2a0afb..ad09d2bf 100644 --- a/tests/pipelines/test_team_98BT.py +++ b/tests/pipelines/test_team_98BT.py @@ -20,8 +20,6 @@ from narps_open.utils.configuration import Configuration from narps_open.pipelines.team_98BT import PipelineTeam98BT -TEMPORARY_DIR = join(Configuration()['directories']['test_runs'], 'test_98BT') - class TestPipelinesTeam98BT: """ A class that contains all the unit tests for the PipelineTeam98BT class.""" @@ -84,8 +82,7 @@ def test_fieldmap_info(): @staticmethod @mark.unit_test - @mark.parametrize('remove_test_dir', TEMPORARY_DIR) - def test_parameters_files(remove_test_dir): + def test_parameters_files(temporary_data_dir): """ Test the get_parameters_files method For this test, we created the two following files by downsampling output files from the preprocessing pipeline : @@ -105,11 +102,11 @@ def test_parameters_files(remove_test_dir): # Get new parameters file PipelineTeam98BT.get_parameters_file( - parameters_file, wc2_file, func_file, 'sid', 'rid', TEMPORARY_DIR) + parameters_file, wc2_file, func_file, 'sid', 'rid', temporary_data_dir) # Check parameters file was created created_parameters_file = join( - TEMPORARY_DIR, 'parameters_files', 'parameters_file_sub-sid_run-rid.tsv') + temporary_data_dir, 'parameters_files', 'parameters_file_sub-sid_run-rid.tsv') assert exists(created_parameters_file) # Check contents From fb023412e0d61f5d226a53e5a7e90bfce8470cda Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= <117362283+bclenet@users.noreply.github.com> Date: Mon, 18 Mar 2024 10:10:45 +0100 Subject: [PATCH 3/8] R9K3 reproduction (#181) * [R9K3][TEST] code and tests for the pipeline * [BUG] smoothing outputs * [BUG] contrast definitions * [BUG] input func from derivs * Select file template naming error * Correcting hypothesis files * Hypotheses files names --- narps_open/pipelines/__init__.py | 2 +- narps_open/pipelines/team_R9K3.py | 638 ++++++++++++++++ narps_open/pipelines/team_R9K3_wip.py | 715 ------------------ tests/pipelines/test_team_R9K3.py | 123 +++ .../pipelines/team_R9K3/confounds.tsv | 3 + 5 files changed, 765 insertions(+), 716 deletions(-) create mode 100644 narps_open/pipelines/team_R9K3.py delete mode 100755 narps_open/pipelines/team_R9K3_wip.py create mode 100644 tests/pipelines/test_team_R9K3.py create mode 100644 tests/test_data/pipelines/team_R9K3/confounds.tsv diff --git a/narps_open/pipelines/__init__.py b/narps_open/pipelines/__init__.py index a9b62d9b..b8ca93d4 100644 --- a/narps_open/pipelines/__init__.py +++ b/narps_open/pipelines/__init__.py @@ -66,7 +66,7 @@ 'R42Q': None, 'R5K7': None, 'R7D1': None, - 'R9K3': None, + 'R9K3': 'PipelineTeamR9K3', 'SM54': None, 'T54A': 'PipelineTeamT54A', 'U26C': 'PipelineTeamU26C', diff --git a/narps_open/pipelines/team_R9K3.py b/narps_open/pipelines/team_R9K3.py new file mode 100644 index 00000000..2946094b --- /dev/null +++ b/narps_open/pipelines/team_R9K3.py @@ -0,0 +1,638 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Write the work of NARPS' team R9K3 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.utility.base import Merge +from nipype.interfaces.io import SelectFiles, DataSink +from nipype.interfaces.spm import ( + Smooth, Level1Design, OneSampleTTestDesign, TwoSampleTTestDesign, + EstimateModel, EstimateContrast, 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 PipelineTeamR9K3(Pipeline): + """ A class that defines the pipeline of team R9K3. """ + + def __init__(self): + super().__init__() + self.fwhm = 6.0 + self.team_id = 'R9K3' + self.contrast_list = ['0001', '0002'] + conditions = ['trialxgain^1', 'trialxloss^1'] + self.subject_level_contrasts = [ + ['effect_of_gain', 'T', conditions, [1, 0]], + ['effect_of_loss', 'T', conditions, [0, 1]] + ] + + def get_preprocessing(self): + """ + Create the preprocessing workflow. + + Returns: + - preprocessing : nipype.WorkFlow + """ + # Initialize preprocessing workflow to connect nodes along the way + preprocessing = Workflow( + base_dir = self.directories.working_dir, + name = 'preprocessing') + + # IDENTITY INTERFACE - To iterate on subjects + information_source = Node(IdentityInterface( + fields = ['subject_id']), + name = 'information_source') + information_source.iterables = [('subject_id', self.subject_list)] + + # SELECT FILES - to select necessary files + templates = { + 'func' : join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-*_bold_space-MNI152NLin2009cAsym_preproc.nii.gz') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory = self.directories.dataset_dir + preprocessing.connect(information_source, 'subject_id', select_files, 'subject_id') + + # GUNZIP - gunzip files because SPM do not use .nii.gz files + gunzip_func = MapNode(Gunzip(), name = 'gunzip_func', iterfield = ['in_file']) + preprocessing.connect(select_files, 'func', gunzip_func, 'in_file') + + # SMOOTH - Spatial smoothing of fMRI data. + smoothing = MapNode(Smooth(), name = 'smoothing', iterfield = 'in_files') + smoothing.inputs.fwhm = [self.fwhm] * 3 + preprocessing.connect(gunzip_func, 'out_file', smoothing, 'in_files') + + # DATA SINK - store the wanted results in the wanted repository + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + preprocessing.connect(smoothing, 'smoothed_files', data_sink, 'preprocessing.@smoothed') + + # Remove large files, if requested + if Configuration()['pipelines']['remove_unused_data']: + + # MERGE - Merge all temporary outputs once they are no longer needed + merge_temp_files = Node(Merge(2), name = 'merge_temp_files') + preprocessing.connect(gunzip_func, 'out_file', merge_temp_files, 'in1') + preprocessing.connect(smoothing, 'smoothed_files', merge_temp_files, 'in2') + + # FUNCTION - Remove temporary files once they are no longer needed + remove_gunziped = MapNode( + InterfaceFactory.create('remove_parent_directory'), + name = 'remove_gunziped', + iterfield = 'file_name' + ) + preprocessing.connect(merge_temp_files, 'out', remove_gunziped, 'file_name') + preprocessing.connect(data_sink, 'out_file', remove_gunziped, '_') + + return preprocessing + + def get_preprocessing_outputs(self): + """ Return the names of the files the preprocessing is supposed to generate. """ + + output_dir = join(self.directories.output_dir, 'preprocessing', '_subject_id_{subject_id}') + + # Smoothing outputs + templates = [join(output_dir, f'_smoothing{index}', + 'ssub-{subject_id}'+f'_task-MGT_run-{run_id}_bold.nii')\ + for index, run_id in zip(range(len(self.run_list)), self.run_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_run_level_analysis(self): + """ No run level analysis has been done by team R9K3 """ + return None + + def get_subject_information(event_file): + """ + Create Bunchs for specifySPMModel, from data extracted from an event_file. + + Parameters : + - event_files: str, event file (one per run) for the subject + + Returns : + - subject_information: Bunch, relevant event information for subject level analysis. + """ + from nipype.interfaces.base import Bunch + from numpy import max as npmax + + trial_onsets = [] + trial_durations = [] + weights_gain = [] + weights_loss = [] + + with open(event_file, 'rt') as file: + next(file) # skip the header + + for line in file: + info = line.strip().split() + + trial_onsets.append(float(info[0])) + trial_durations.append(0.0) + weights_gain.append(float(info[2])) + weights_loss.append(float(info[3])) + + # Scale weights between 0 and 1 + weights_gain = list(weights_gain / npmax(weights_gain)) + weights_loss = list(weights_loss / npmax(weights_loss)) + + return Bunch( + conditions = ['trial'], + onsets = [trial_onsets], + durations = [trial_durations], + amplitudes = None, + tmod = None, + pmod = [ + Bunch( + name = ['gain', 'loss'], + poly = [1, 1], + param = [weights_gain, weights_loss], + ), + None, + ], + regressor_names = None, + regressors = None + ) + + def get_confounds_file(confounds_file: str, subject_id: str, run_id: str) -> str: + """ + Create a tsv file with only desired confounds per subject per run. + + Parameters : + - confounds_file: str, path to the file containing confounds from fmriprep + - subject_id : related subject id + - run_id : related run id + + Return : + - out_file : path to new file containing only desired confounds + """ + from os.path import abspath + + from pandas import read_csv, DataFrame + from numpy import array, transpose + + # Get the dataframe containing the 6 head motion parameter regressors + data_frame = read_csv(confounds_file, sep = '\t', header=0) + + # Extract parameters 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']]))) + + # Write confounds to a file + out_file = abspath(f'confounds_file_sub-{subject_id}_run-{run_id}.tsv') + with open(out_file, 'w', encoding = 'utf-8') as writer: + writer.write(retained_parameters.to_csv( + sep = '\t', index = False, header = False, na_rep = '0.0')) + + return out_file + + def get_subject_level_analysis(self): + """ + Create the subject level analysis workflow. + + Returns: + - subject_level : nipype.WorkFlow + WARNING: the name attribute of the workflow is 'subject_level_analysis' + """ + # Create subject level analysis workflow + subject_level = Workflow( + base_dir = self.directories.working_dir, + name = 'subject_level_analysis') + + # IDENTITY INTERFACE - To iterate on subjects + information_source = Node(IdentityInterface( + fields = ['subject_id']), + name = 'information_source') + information_source.iterables = [('subject_id', self.subject_list)] + + # SELECT FILES - to select necessary files + templates = { + 'confounds' : join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-*_bold_confounds.tsv'), + 'func' : join(self.directories.output_dir, 'preprocessing', '_subject_id_{subject_id}', + '_smoothing*', + 'ssub-{subject_id}_task-MGT_run-*_bold_space-MNI152NLin2009cAsym_preproc.nii'), + 'event' : 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 + subject_level.connect(information_source, 'subject_id', select_files, 'subject_id') + + # FUNCTION get_subject_information - generate files with event data + subject_information = MapNode(Function( + function = self.get_subject_information, + input_names = ['event_file'], + output_names = ['subject_info']), + iterfield = 'event_file', + name = 'subject_information') + subject_level.connect(select_files, 'event', subject_information, 'event_file') + + # FUNCTION node get_confounds_file - generate files with confounds data + confounds = MapNode( + Function( + function = self.get_confounds_file, + input_names = ['confounds_file', 'subject_id', 'run_id'], + output_names = ['confounds_file'] + ), + name = 'confounds', + iterfield = ['confounds_file', 'run_id']) + confounds.inputs.run_id = self.run_list + subject_level.connect(information_source, 'subject_id', confounds, 'subject_id') + subject_level.connect(select_files, 'confounds', confounds, 'confounds_file') + + # SPECIFY MODEL - generates SPM-specific Model + specify_model = Node(SpecifySPMModel(), name = 'specify_model') + 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 + subject_level.connect(select_files, 'func', specify_model, 'functional_runs') + subject_level.connect(confounds, 'confounds_file', specify_model, 'realignment_parameters') + subject_level.connect(subject_information, 'subject_info', specify_model, 'subject_info') + + # LEVEL 1 DESIGN - Generates an SPM design matrix + model_design = Node(Level1Design(), name = 'model_design') + model_design.inputs.bases = {'hrf': {'derivs': [0, 0]}} + model_design.inputs.timing_units = 'secs' + model_design.inputs.interscan_interval = TaskInformation()['RepetitionTime'] + subject_level.connect(specify_model, 'session_info', model_design, 'session_info') + + # ESTIMATE MODEL - estimate the parameters of the model + model_estimate = Node(EstimateModel(), name = 'model_estimate') + model_estimate.inputs.estimation_method = {'Classical': 1} + subject_level.connect(model_design, 'spm_mat_file', model_estimate, 'spm_mat_file') + + # ESTIMATE CONTRAST - estimates contrasts + contrast_estimate = Node(EstimateContrast(), name = 'contrast_estimate') + contrast_estimate.inputs.contrasts = self.subject_level_contrasts + subject_level.connect([ + (model_estimate, contrast_estimate, [ + ('spm_mat_file', 'spm_mat_file'), + ('beta_images', 'beta_images'), + ('residual_image', 'residual_image') + ]) + ]) + + # DATA SINK - 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, data_sink, [ + ('con_images', 'subject_level_analysis.@con_images'), + ('spmT_images', 'subject_level_analysis.@spmT_images'), + ('spm_mat_file', 'subject_level_analysis.@spm_mat_file') + ]) + ]) + + return subject_level + + def get_subject_level_outputs(self): + """ Return the names of the files the subject level analysis is supposed to generate. """ + + # Contrat maps + templates = [join( + self.directories.output_dir, + 'subject_level_analysis', '_subject_id_{subject_id}', f'con_{contrast_id}.nii')\ + for contrast_id in self.contrast_list] + + # SPM.mat file + templates += [join( + self.directories.output_dir, + 'subject_level_analysis', '_subject_id_{subject_id}', 'SPM.mat')] + + # spmT maps + templates += [join( + self.directories.output_dir, + 'subject_level_analysis', '_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 + """ + + 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_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 + information_source = Node( + IdentityInterface( + fields=['contrast_id']), + name='information_source') + information_source.iterables = [('contrast_id', self.contrast_list)] + + # SelectFiles + templates = { + # Contrasts for all participants + 'contrasts' : join(self.directories.output_dir, + 'subject_level_analysis', '_subject_id_*', 'con_{contrast_id}.nii') + } + + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory = self.directories.results_dir + select_files.inputs.force_lists = True + + # Datasink - save important files + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + + # 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 + + # 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' + ) + + # Estimate model + estimate_model = Node(EstimateModel(), name = 'estimate_model') + estimate_model.inputs.estimation_method = {'Classical':1} + + # Estimate contrasts + estimate_contrast = Node(EstimateContrast(), name = 'estimate_contrast') + estimate_contrast.inputs.group_contrast = True + + ## Create thresholded maps + threshold = MapNode(Threshold(), + name = 'threshold', + iterfield = ['stat_image', 'contrast_index']) + threshold.inputs.use_fwe_correction = False + threshold.inputs.height_threshold = 0.001 + threshold.inputs.extent_threshold = 5 + threshold.synchronize = True + + group_level_analysis = Workflow( + base_dir = self.directories.working_dir, + name = f'group_level_analysis_{method}_nsub_{nb_subjects}') + group_level_analysis.connect([ + (information_source, select_files, [('contrast_id', 'contrast_id')]), + (select_files, get_contrasts, [('contrasts', 'input_str')]), + (estimate_model, estimate_contrast, [ + ('spm_mat_file', 'spm_mat_file'), + ('residual_image', 'residual_image'), + ('beta_images', 'beta_images')]), + (estimate_contrast, threshold, [ + ('spm_mat_file', 'spm_mat_file'), + ('spmT_images', 'stat_image')]), + (estimate_model, data_sink, [ + ('mask_image', f'group_level_analysis_{method}_nsub_{nb_subjects}.@mask')]), + (estimate_contrast, data_sink, [ + ('spm_mat_file', f'group_level_analysis_{method}_nsub_{nb_subjects}.@spm_mat'), + ('spmT_images', f'group_level_analysis_{method}_nsub_{nb_subjects}.@T'), + ('con_images', f'group_level_analysis_{method}_nsub_{nb_subjects}.@con')]), + (threshold, data_sink, [ + ('thresholded_map', f'group_level_analysis_{method}_nsub_{nb_subjects}.@thresh')])]) + + if method in ('equalRange', 'equalIndifference'): + estimate_contrast.inputs.contrasts = [ + ('Group', 'T', ['mean'], [1]), + ('Group', 'T', ['mean'], [-1]) + ] + + threshold.inputs.contrast_index = [1, 2] + + # Specify design matrix + one_sample_t_test_design = Node(OneSampleTTestDesign(), + name = 'one_sample_t_test_design') + group_level_analysis.connect([ + (get_contrasts, one_sample_t_test_design, [ + (('out_list', clean_list), 'in_files') + ]), + (one_sample_t_test_design, estimate_model, [('spm_mat_file', 'spm_mat_file')]) + ]) + + if method == 'equalRange': + group_level_analysis.connect([ + (get_equal_range_subjects, get_contrasts, [ + (('out_list', complete_subject_ids), 'elements') + ]) + ]) + + elif method == 'equalIndifference': + group_level_analysis.connect([ + (get_equal_indifference_subjects, get_contrasts, [ + (('out_list', complete_subject_ids), 'elements') + ]) + ]) + + elif method == 'groupComp': + estimate_contrast.inputs.contrasts = [ + ('Eq range vs Eq indiff in loss', 'T', ['Group_{1}', 'Group_{2}'], [1, -1]) + ] + + threshold.inputs.contrast_index = [1] + + # 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_2 = MapNode(Function( + function = elements_in_string, + input_names = ['input_str', 'elements'], + output_names = ['out_list'] + ), + name = 'get_contrasts_2', iterfield = 'input_str' + ) + + # Specify design matrix + two_sample_t_test_design = Node(TwoSampleTTestDesign(), + name = 'two_sample_t_test_design') + two_sample_t_test_design.inputs.unequal_variance = True + + group_level_analysis.connect([ + (select_files, get_contrasts_2, [('contrasts', 'input_str')]), + (get_equal_range_subjects, get_contrasts, [ + (('out_list', complete_subject_ids), 'elements') + ]), + (get_equal_indifference_subjects, get_contrasts_2, [ + (('out_list', complete_subject_ids), 'elements') + ]), + (get_contrasts, two_sample_t_test_design, [ + (('out_list', clean_list), 'group1_files') + ]), + (get_contrasts_2, two_sample_t_test_design, [ + (('out_list', clean_list), 'group2_files') + ]), + (two_sample_t_test_design, estimate_model, [('spm_mat_file', 'spm_mat_file')]) + ]) + + 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. + Note that hypotheses 5 to 8 correspond to the maps given by the team in their results ; + but they are not fully consistent with the hypotheses definitions as expected by NARPS. + """ + 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/narps_open/pipelines/team_R9K3_wip.py b/narps_open/pipelines/team_R9K3_wip.py deleted file mode 100755 index 65613e70..00000000 --- a/narps_open/pipelines/team_R9K3_wip.py +++ /dev/null @@ -1,715 +0,0 @@ -# pylint: skip-file -# THIS IS A TEMPLATE THAT CAN BE USE TO REPRODUCE A NEW PIPELINE - -import os -import shutil -from os.path import join as opj -from typing import List - -from nipype import Node, Workflow -from nipype.algorithms.misc import Gunzip -from nipype.algorithms.modelgen import ( # Functions used during L1 analysis - SpecifyModel, - SpecifySPMModel, -) -from nipype.interfaces.base import Bunch -from nipype.interfaces.io import DataSink, SelectFiles -from nipype.interfaces.spm import Smooth -from nipype.interfaces.utility import Function, IdentityInterface - -from .utils import fmriprep_data_template, raw_data_template - - -def get_preprocessing( - exp_dir: str, - result_dir: str, - working_dir: str, - output_dir: str, - subject_list: List[str], - run_list: List[str], - fwhm: float, -): - """ - Returns the preprocessing workflow. - - Parameters: - - exp_dir: directory where raw data are stored - - result_dir: directory where results will be stored - - working_dir: name of the sub-directory for intermediate results - - output_dir: name of the sub-directory for final results - - subject_list: list of subject for which you want to do the preprocessing - - run_list: list of runs for which you want to do the preprocessing - - fwhm: fwhm for smoothing step - - Returns: - - preprocessing: Nipype WorkFlow - """ - - infosource_preproc = Node( - IdentityInterface(fields=["subject_id", "run_id"]), name="infosource_preproc" - ) - - # Iterates over subject and runs - infosource_preproc.iterables = [ - ("subject_id", subject_list), - ("run_id", run_list), - ] - - # SelectFiles node - to select necessary files - selectfiles_preproc = Node( - SelectFiles(fmriprep_data_template(), base_directory=exp_dir), - name="selectfiles_preproc", - ) - - # DataSink Node - store the wanted results in the wanted repository - datasink_preproc = Node( - DataSink(base_directory=result_dir, container=working_dir), - name="datasink_preproc", - ) - - gunzip_func = Node(Gunzip(), name="gunzip_func") - - smooth = Node(Smooth(fwhm=fwhm, implicit_masking=False), name="smooth") - - preprocessing = Workflow( - base_dir=opj(result_dir, working_dir), name="preprocessing" - ) - - preprocessing.connect( - [ - ( - infosource_preproc, - selectfiles_preproc, - [("subject_id", "subject_id"), ("run_id", "run_id")], - ), - ( - selectfiles_preproc, - gunzip_func, - [("func_preproc", "in_file")], - ), - ( - gunzip_func, - smooth, - [("out_file", "in_files")], - ), - ( - smooth, - datasink_preproc, - [("smoothed_files", "preprocess.@sym_link")], - ), - ] - ) - - return preprocessing - -""" -In first level analyses, three regressors of interest were included in the model: - -1) task-related activity, -2) task-related activity modulated by the amount of gain, and -3) task-related activity modulated by the amount of loss. - -The regressor for the basic task-related activity was defined -as a vector of ones at every trial onset and zeros at all other time points -(event duration = 0 s). - -Within each functional run, each parametric modulator -(the raw gain or loss values) was scaled such that its maximum value became 1. -All three regressors were convolved with the canonical hemodynamic response function. -The two modulated regressors were orthogonalized against the original task regressor. -Temporal/dispersion derivatives of the regressors were not included. - -In both first and second level analyses, model parameters were estimated -using the classical (Restricted Maximum Likelihood) method as implemented in SPM12. - -For second level analyses, we used a mixed-effects approach with weighted least squares. - -In the second level analysis performed to test hypothesis 9 -(two-sample t-test comparing the equal indifference group and the equal range group), -the two groups were assumed to be independent and have unequal variance. -""" - - -# FIXME: THIS FUNCTION IS USED IN THE FIRST LEVEL ANALYSIS PIPELINES OF SPM -# THIS IS AN EXAMPLE THAT IS ADAPTED TO A SPECIFIC PIPELINE -# MODIFY ACCORDING TO THE PIPELINE YOU WANT TO REPRODUCE -def get_subject_infos_spm(event_files: List[str], runs: List[str]): - """ - Create Bunchs for specifySPMModel. - - Parameters : - - event_files: list of events files (one per run) for the subject - - runs: list of runs to use - - Returns : - - subject_info : list of Bunch for 1st level analysis. - """ - - cond_names = ["trial", "accepting", "rejecting"] - onset = {} - duration = {} - weights_gain = {} - weights_loss = {} - onset_button = {} - duration_button = {} - - # Loop over number of runs. - for r in range(len(runs)): - onset |= {f"{s}_run{str(r + 1)}": [] for s in cond_names} - duration |= {f"{s}_run{str(r + 1)}": [] for s in cond_names} - weights_gain[f"gain_run{str(r + 1)}"] = [] - weights_loss[f"loss_run{str(r + 1)}"] = [] - - for r, run in enumerate(runs): - - f_events = event_files[r] - - with open(f_events, "rt") as f: - next(f) # skip the header - - for line in f: - info = line.strip().split() - - for cond in cond_names: - val = f"{cond}_run{str(r + 1)}" - val_gain = f"gain_run{str(r + 1)}" - val_loss = f"loss_run{str(r + 1)}" - if cond == "trial": - onset[val].append(float(info[0])) # onsets for trial_run1 - duration[val].append(float(4)) - weights_gain[val_gain].append( - float(info[2]) - ) # weights gain for trial_run1 - weights_loss[val_loss].append( - float(info[3]) - ) # weights loss for trial_run1 - elif cond == "accepting" and "accept" in info[5]: - onset[val].append(float(info[0]) + float(info[4])) - duration[val].append(float(0)) - elif cond == "rejecting" and "reject" in info[5]: - onset[val].append(float(info[0]) + float(info[4])) - duration[val].append(float(0)) - - # Bunching is done per run, i.e. trial_run1, trial_run2, etc. - # But names must not have '_run1' etc because we concatenate runs - subject_info = [] - for r in range(len(runs)): - - cond = [f"{s}_run{str(r + 1)}" for s in cond_names] - gain = f"gain_run{str(r + 1)}" - loss = f"loss_run{str(r + 1)}" - - subject_info.insert( - r, - Bunch( - conditions=cond_names, - onsets=[onset[c] for c in cond], - durations=[duration[c] for c in cond], - amplitudes=None, - tmod=None, - pmod=[ - Bunch( - name=["gain", "loss"], - poly=[1, 1], - param=[weights_gain[gain], weights_loss[loss]], - ), - None, - ], - regressor_names=None, - regressors=None, - ), - ) - - return subject_info - - -# FIXME: THIS FUNCTION CREATES THE CONTRASTS THAT WILL BE ANALYZED IN THE FIRST LEVEL ANALYSIS -# IT IS ADAPTED FOR A SPECIFIC PIPELINE AND SHOULD BE MODIFIED DEPENDING ON THE PIPELINE -# YOU ARE TRYING TO REPRODUCE -def get_contrasts(subject_id: str): - """ - Create the list of tuples that represents contrasts. - Each contrast is in the form : - (Name,Stat,[list of condition names],[weights on those conditions]) - - Parameters: - - subject_id: ID of the subject - - Returns: - - contrasts: list of tuples, list of contrasts to analyze - """ - # list of condition names - conditions = ["trial", "trialxgain^1", "trialxloss^1"] - - # create contrasts - trial = ("trial", "T", conditions, [1, 0, 0]) - - effect_gain = ("effect_of_gain", "T", conditions, [0, 1, 0]) - - effect_loss = ("effect_of_loss", "T", conditions, [0, 0, 1]) - - # contrast list - contrasts = [effect_gain, effect_loss] - - return contrasts - - -# FUNCTION TO CREATE THE WORKFLOW OF A L1 ANALYSIS (SUBJECT LEVEL) -def get_l1_analysis( - exp_dir: str, - result_dir: str, - working_dir: str, - output_dir: str, - subject_list: List[str], - run_list: List[str], - TR: float, -): - """ - Returns the first level analysis workflow. - - Parameters: - - exp_dir: directory where raw data are stored - - result_dir: directory where results will be stored - - working_dir: name of the sub-directory for intermediate results - - output_dir: name of the sub-directory for final results - - subject_list: list of subject for which you want to do the analysis - - run_list: list of runs for which you want to do the analysis - - TR: time repetition used during acquisition - - Returns: - - l1_analysis : Nipype WorkFlow - """ - # THE FOLLOWING PART STAYS THE SAME FOR ALL PIPELINES - # Infosource Node - To iterate on subjects - infosource = Node( - IdentityInterface( - fields=["subject_id", "exp_dir", "result_dir", "working_dir", "run_list"], - exp_dir=exp_dir, - result_dir=result_dir, - working_dir=working_dir, - run_list=run_list, - ), - name="infosource", - ) - - # ITERATES OVER SUBJECT LIST - infosource.iterables = [("subject_id", subject_list)] - - # Templates to select files node - - # FIXME: CHANGE THE NAME OF THE FILE - # DEPENDING ON THE FILENAMES OF RESULTS OF PREPROCESSING - func_file = opj( - result_dir, - output_dir, - "preprocess", - "_run_id_*_subject_id_{subject_id}", - "complete_filename_{subject_id}_complete_filename.nii", - ) - - event_files = opj( - exp_dir, - "sub-{subject_id}", - "func", - "sub-{subject_id}_task-MGT_run-*_events.tsv", - ) - - - - template = {"func": func_file, "event": raw_data_template()["event_file"]} - - # SelectFiles node - to select necessary files - selectfiles = Node( - SelectFiles(template, base_directory=exp_dir), name="selectfiles" - ) - - # DataSink Node - store the wanted results in the wanted repository - datasink = Node( - DataSink(base_directory=result_dir, container=output_dir), name="datasink" - ) - - # Get Subject Info - get subject specific condition information - subject_infos = Node( - Function( - input_names=["event_files", "runs"], - output_names=["subject_info"], - function=get_subject_infos_spm, - ), - name="subject_infos", - ) - - subject_infos.inputs.runs = run_list - # THIS IS THE NODE EXECUTING THE get_contrasts FUNCTION - # Node contrasts to get contrasts - contrasts = Node( - Function( - function=get_contrasts, - input_names=["subject_id"], - output_names=["contrasts"], - ), - name="contrasts", - ) - - # Create l1 analysis workflow and connect its nodes - l1_analysis = Workflow(base_dir=opj(result_dir, working_dir), name="l1_analysis") - - l1_analysis.connect( - [ - (infosource, selectfiles, [("subject_id", "subject_id")]), - (infosource, contrasts, [("subject_id", "subject_id")]), - (selectfiles, subject_infos, [("event", "event_files")]), - # FIXME: Complete with name of node to link with and the name of the input - ( - selectfiles, - node_variable[("func", "node_input_name")], - ), - # Input and output names can be found on NiPype documentation - (node_variable, datasink, [("node_output_name", "preprocess.@sym_link")]), - ] - ) - - return l1_analysis - - -# THIS FUNCTION RETURNS THE LIST OF IDS AND FILES OF EACH GROUP OF PARTICIPANTS -# TO DO SEPARATE GROUP LEVEL ANALYSIS AND BETWEEN GROUP ANALYSIS -# THIS FUNCTIONS IS ADAPTED FOR AN SPM PIPELINE. -def get_subset_contrasts_spm( - file_list, subject_list: List[str], participants_file: str -): - """ - Parameters : - - file_list : original file list selected by selectfiles node - - subject_list : list of subject IDs that are in the wanted group for the analysis - - participants_file: file containing participants characteristics - - This function return the file list containing only the files belonging - to the subject in the wanted group. - """ - equalIndifference_id = [] - equalRange_id = [] - equalIndifference_files = [] - equalRange_files = [] - - with open( - participants_file, "rt" - ) as f: # Reading file containing participants IDs and groups - next(f) # skip the header - - for line in f: - info = line.strip().split() - - if info[0][-3:] in subject_list: - if info[1] == "equalIndifference": # Checking for each participant if its ID was selected - # and separate people depending on their group - equalIndifference_id.append(info[0][-3:]) - elif info[1] == "equalRange": - equalRange_id.append(info[0][-3:]) - - # Checking for each selected file if the corresponding participant was selected - # and add the file to the list corresponding to its group - for file in file_list: - sub_id = file.split("/") - if sub_id[-2][-3:] in equalIndifference_id: - equalIndifference_files.append(file) - elif sub_id[-2][-3:] in equalRange_id: - equalRange_files.append(file) - - return ( - equalIndifference_id, - equalRange_id, - equalIndifference_files, - equalRange_files, - ) - - -# FUNCTION TO CREATE THE WORKFLOW OF A L2 ANALYSIS (GROUP LEVEL) -def get_l2_analysis( - exp_dir: str, - result_dir: str, - working_dir: str, - output_dir: str, - subject_list: List[str], - contrast_list: List[str], - n_sub: int, - method: str, -): - """ - Returns the 2nd level of analysis workflow. - - Parameters: - - exp_dir: directory where raw data are stored - - result_dir: directory where results will be stored - - working_dir: name of the sub-directory for intermediate results - - output_dir: name of the sub-directory for final results - - subject_list: list of subject for which you want to do the preprocessing - - contrast_list: list of contrasts to analyze - - n_sub: number of subjects used to do the analysis - - method: one of "equalRange", "equalIndifference" or "groupComp" - - Returns: - - l2_analysis: Nipype WorkFlow - """ - # THE FOLLOWING PART STAYS THE SAME FOR ALL PREPROCESSING PIPELINES - # Infosource - a function free node to iterate over the list of subject names - infosource_groupanalysis = Node( - IdentityInterface(fields=["contrast_id", "subjects"], subjects=subject_list), - name="infosource_groupanalysis", - ) - - infosource_groupanalysis.iterables = [("contrast_id", contrast_list)] - - # SelectFiles - contrast_file = opj( - result_dir, - output_dir, - "l1_analysis", - "_subject_id_*", - "complete_filename_{contrast_id}_complete_filename.nii", - ) - # FIXME: CHANGE THE NAME OF THE FILE DEPENDING ON - # THE FILENAMES OF THE RESULTS OF PREPROCESSING - # (DIFFERENT FOR AN FSL PIPELINE) - - participants_file = opj(exp_dir, "participants.tsv") - - templates = {"contrast": contrast_file, "participants": participants_file} - - selectfiles_groupanalysis = Node( - SelectFiles(templates, base_directory=result_dir, force_list=True), - name="selectfiles_groupanalysis", - ) - - # Datasink node : to save important files - datasink_groupanalysis = Node( - DataSink(base_directory=result_dir, container=output_dir), - name="datasink_groupanalysis", - ) - - # IF THIS IS AN SPM PIPELINE: - # Node to select subset of contrasts - sub_contrasts = Node( - Function( - input_names=["file_list", "method", "subject_list", "participants_file"], - output_names=[ - "equalIndifference_id", - "equalRange_id", - "equalIndifference_files", - "equalRange_files", - ], - function=get_subset_contrasts_spm, - ), - name="sub_contrasts", - ) - - sub_contrasts.inputs.method = method - - # IF THIS IS AN FSL PIPELINE: - subgroups_contrasts = Node( - Function( - input_names=["copes", "varcopes", "subject_ids", "participants_file"], - output_names=[ - "copes_equalIndifference", - "copes_equalRange", - "varcopes_equalIndifference", - "varcopes_equalRange", - "equalIndifference_id", - "equalRange_id", - "copes_global", - "varcopes_global", - ], - function=get_subgroups_contrasts, - ), - name="subgroups_contrasts", - ) - - regs = Node( - Function( - input_names=[ - "equalRange_id", - "equalIndifference_id", - "method", - "subject_list", - ], - output_names=["regressors"], - function=get_regs, - ), - name="regs", - ) - regs.inputs.method = method - regs.inputs.subject_list = subject_list - - # FIXME: THE FOLLOWING PART HAS TO BE MODIFIED WITH NODES OF THE PIPELINE - node_variable = Node( - node_function, name="node_name" - ) # Replace with the name of the node_variable, - # the node_function to use in the NiPype interface, - # and the name of the node (recommended to be the same as node_variable) - - # FIXME: ADD OTHER NODES WITH THE DIFFERENT STEPS OF THE PIPELINE - - l2_analysis = Workflow( - base_dir=opj(result_dir, working_dir), name=f"l2_analysis_{method}_nsub_{n_sub}" - ) - # FOR AN SPM PIPELINE - l2_analysis.connect( - [ - ( - infosource_groupanalysis, - selectfiles_groupanalysis, - [("contrast_id", "contrast_id")], - ), - (infosource_groupanalysis, sub_contrasts, [("subjects", "subject_list")]), - ( - selectfiles_groupanalysis, - sub_contrasts, - [("contrast", "file_list"), ("participants", "participants_file")], - ), # Complete with other links between nodes - ] - ) - - # FOR AN FSL PIPELINE - l2_analysis.connect( - [ - ( - infosource_groupanalysis, - selectfiles_groupanalysis, - [("contrast_id", "contrast_id")], - ), - ( - infosource_groupanalysis, - subgroups_contrasts, - [("subject_list", "subject_ids")], - ), - ( - selectfiles_groupanalysis, - subgroups_contrasts, - [ - ("cope", "copes"), - ("varcope", "varcopes"), - ("participants", "participants_file"), - ], - ), - ( - selectfiles_groupanalysis, - node_variable[("func", "node_input_name")], - ), # Complete with name of node to link with and the name of the input - # Input and output names can be found on NiPype documentation - ( - node_variable, - datasink_groupanalysis, - [("node_output_name", "preprocess.@sym_link")], - ), - ] - ) # Complete with other links between nodes - - if method == "equalRange" or method == "equalIndifference": - contrasts = [("Group", "T", ["mean"], [1]), ("Group", "T", ["mean"], [-1])] - - elif method == "groupComp": - contrasts = [ - ("Eq range vs Eq indiff in loss", "T", ["Group_{1}", "Group_{2}"], [1, -1]) - ] - - # FIXME: ADD OTHER NODES WITH THE DIFFERENT STEPS OF THE PIPELINE - - return l2_analysis - - -# THIS FUNCTION IS USED TO REORGANIZE FINAL RESULTS OF THE PIPELINE -def reorganize_results(result_dir: str, output_dir: str, n_sub: int, team_ID: str): - """ - Reorganize the results to analyze them. - - Parameters: - - result_dir: directory where results will be stored - - output_dir: name of the sub-directory for final results - - n_sub: number of subject used for the analysis - - team_ID: ID of the team to reorganize results - - """ - - h1 = opj( - result_dir, - output_dir, - f"l2_analysis_equalIndifference_nsub_{n_sub}", - "_contrast_id_01", - ) - h2 = opj( - result_dir, - output_dir, - f"l2_analysis_equalRange_nsub_{n_sub}", - "_contrast_id_01", - ) - h3 = opj( - result_dir, - output_dir, - f"l2_analysis_equalIndifference_nsub_{n_sub}", - "_contrast_id_01", - ) - h4 = opj( - result_dir, - output_dir, - f"l2_analysis_equalRange_nsub_{n_sub}", - "_contrast_id_01", - ) - h5 = opj( - result_dir, - output_dir, - f"l2_analysis_equalIndifference_nsub_{n_sub}", - "_contrast_id_02", - ) - h6 = opj( - result_dir, - output_dir, - f"l2_analysis_equalRange_nsub_{n_sub}", - "_contrast_id_02", - ) - h7 = opj( - result_dir, - output_dir, - f"l2_analysis_equalIndifference_nsub_{n_sub}", - "_contrast_id_02", - ) - h8 = opj( - result_dir, - output_dir, - f"l2_analysis_equalRange_nsub_{n_sub}", - "_contrast_id_02", - ) - h9 = opj( - result_dir, output_dir, f"l2_analysis_groupComp_nsub_{n_sub}", "_contrast_id_02" - ) - - h = [h1, h2, h3, h4, h5, h6, h7, h8, h9] - - repro_unthresh = [ - opj(filename, "_change_filename_.nii") for i, filename in enumerate(h) - ] # Change filename with the filename of the final results - - repro_thresh = [ - opj(filename, "_change_filename_.nii") for i, filename in enumerate(h) - ] - - if not os.path.isdir(opj(result_dir, "NARPS-reproduction")): - os.mkdir(opj(result_dir, "NARPS-reproduction")) - - for i, filename in enumerate(repro_unthresh): - f_in = filename - f_out = opj( - result_dir, - "NARPS-reproduction", - f"team_{team_ID}_nsub_{n_sub}_hypo{i+1}_unthresholded.nii", - ) - shutil.copyfile(f_in, f_out) - - for i, filename in enumerate(repro_thresh): - f_in = filename - f_out = opj( - result_dir, - "NARPS-reproduction", - f"team_{team_ID}_nsub_{n_sub}_hypo{i+1}_thresholded.nii", - ) - shutil.copyfile(f_in, f_out) - - print(f"Results files of team {team_ID} reorganized.") diff --git a/tests/pipelines/test_team_R9K3.py b/tests/pipelines/test_team_R9K3.py new file mode 100644 index 00000000..2a435faf --- /dev/null +++ b/tests/pipelines/test_team_R9K3.py @@ -0,0 +1,123 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Tests of the 'narps_open.pipelines.team_R9K3' module. + +Launch this test with PyTest + +Usage: +====== + pytest -q test_team_R9K3.py + pytest -q test_team_R9K3.py -k +""" +from os.path import join, exists, abspath +from filecmp import cmp + +from pytest import helpers, mark +from nipype import Workflow, Node +from nipype.interfaces.utility import Function +from nipype.interfaces.base import Bunch + +from narps_open.utils.configuration import Configuration +from narps_open.pipelines.team_R9K3 import PipelineTeamR9K3 + +class TestPipelinesTeamR9K3: + """ A class that contains all the unit tests for the PipelineTeamR9K3 class.""" + + @staticmethod + @mark.unit_test + def test_create(): + """ Test the creation of a PipelineTeamR9K3 object """ + + pipeline = PipelineTeamR9K3() + + # 1 - check the parameters + assert pipeline.fwhm == 6.0 + assert pipeline.team_id == 'R9K3' + + # 2 - check workflows + assert isinstance(pipeline.get_preprocessing(), Workflow) + 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 PipelineTeamR9K3 object """ + pipeline = PipelineTeamR9K3() + # 1 - 1 subject outputs + pipeline.subject_list = ['001'] + helpers.test_pipeline_outputs(pipeline, [4, 0, 5, 42, 18]) + + # 2 - 4 subjects outputs + pipeline.subject_list = ['001', '002', '003', '004'] + helpers.test_pipeline_outputs(pipeline, [16, 0, 20, 42, 18]) + + @staticmethod + @mark.unit_test + def test_confounds_file(temporary_data_dir): + """ Test the get_confounds_file method """ + + # Test files + in_confounds_file = abspath(join(Configuration()['directories']['test_data'], + 'pipelines', 'confounds.tsv')) + + # Reference file + ref_file = abspath(join(Configuration()['directories']['test_data'], + 'pipelines', 'team_R9K3', 'confounds.tsv')) + + # Create average values file + confounds_node = Node(Function( + input_names = ['confounds_file', 'subject_id', 'run_id'], + output_names = ['out_file'], + function = PipelineTeamR9K3.get_confounds_file), + name = 'confounds_node') + confounds_node.base_dir = temporary_data_dir + confounds_node.inputs.confounds_file = in_confounds_file + confounds_node.inputs.subject_id = 'sid' + confounds_node.inputs.run_id = 'rid' + confounds_node.run() + + # Check 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) + + # Check contents + assert cmp(ref_file, created_confounds_file) + + @staticmethod + @mark.unit_test + def test_subject_information(): + """ Test the get_subject_information method """ + + # Test with 'gain' + test_event_file = join( + Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') + information = PipelineTeamR9K3.get_subject_information(test_event_file) + + assert isinstance(information, Bunch) + assert information.conditions == ['trial'] + + helpers.compare_float_2d_arrays([[0.0, 0.0, 0.0, 0.0, 0.0]], information.durations) + helpers.compare_float_2d_arrays( + [[4.071, 11.834, 19.535, 27.535, 36.435]], + information.onsets) + paramateric_modulation = information.pmod[0] + assert isinstance(paramateric_modulation, Bunch) + assert paramateric_modulation.name == ['gain', 'loss'] + assert paramateric_modulation.poly == [1, 1] + helpers.compare_float_2d_arrays( + [[0.368421053, 0.894736842, 1.0, 0.263157895, 0.421052632], + [0.315789474, 0.736842105, 1.0, 0.789473684, 0.894736842]], + paramateric_modulation.param) + + @staticmethod + @mark.pipeline_test + def test_execution(): + """ Test the execution of a PipelineTeamR9K3 and compare results """ + helpers.test_pipeline_evaluation('R9K3') diff --git a/tests/test_data/pipelines/team_R9K3/confounds.tsv b/tests/test_data/pipelines/team_R9K3/confounds.tsv new file mode 100644 index 00000000..cf63c178 --- /dev/null +++ b/tests/test_data/pipelines/team_R9K3/confounds.tsv @@ -0,0 +1,3 @@ +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 +-2.56954e-05 -0.00923735 0.0549667 0.000997278 -0.00019745 -0.000398988 From b7f30a18d5a829e4b12aef93bc59b50279b48502 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= <117362283+bclenet@users.noreply.github.com> Date: Thu, 11 Apr 2024 16:59:16 +0200 Subject: [PATCH 4/8] 3TR7 reproduction (#186) * [REPRO][TEST] starting 3TR7 * Issue with gunzip * Regressors naming * Regressors naming * Contrasts naming * Contrasts naming * Hypotheses naming * Test data * Test correaction * Output file names * Output file names --- narps_open/pipelines/__init__.py | 2 +- narps_open/pipelines/team_3TR7.py | 663 ++++++++++++++++++ tests/pipelines/test_team_3TR7.py | 122 ++++ .../pipelines/team_3TR7/confounds.tsv | 3 + 4 files changed, 789 insertions(+), 1 deletion(-) create mode 100644 narps_open/pipelines/team_3TR7.py create mode 100644 tests/pipelines/test_team_3TR7.py create mode 100644 tests/test_data/pipelines/team_3TR7/confounds.tsv diff --git a/narps_open/pipelines/__init__.py b/narps_open/pipelines/__init__.py index b8ca93d4..4294853c 100644 --- a/narps_open/pipelines/__init__.py +++ b/narps_open/pipelines/__init__.py @@ -23,7 +23,7 @@ '2T7P': None, '3C6G': None, '3PQ2': None, - '3TR7': None, + '3TR7': 'PipelineTeam3TR7', '43FJ': None, '46CD': None, '4SZ2': None, diff --git a/narps_open/pipelines/team_3TR7.py b/narps_open/pipelines/team_3TR7.py new file mode 100644 index 00000000..d7d4243d --- /dev/null +++ b/narps_open/pipelines/team_3TR7.py @@ -0,0 +1,663 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Write the work of NARPS' team 3TR7 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 PipelineTeam3TR7(Pipeline): + """ A class that defines the pipeline of team 3TR7. """ + + def __init__(self): + super().__init__() + self.fwhm = 8.0 + self.team_id = '3TR7' + self.contrast_list = ['0001', '0002'] + self.subject_level_contrasts = [ + ('effect_of_gain', 'T', ['gamble', 'gamblexgain^1', 'gamblexloss^1'], [0, 1, 0]), + ('effect_of_loss', 'T', ['gamble', 'gamblexgain^1', 'gamblexloss^1'], [0, 0, 1]) + ] + + def get_preprocessing(self): + """ No preprocessing has been done by team 3TR7 """ + return None + + def get_run_level_analysis(self): + """ No run level analysis has been done by team 3TR7 """ + 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 + + onsets = {} + durations = {} + weights_gain = {} + weights_loss = {} + + subject_info = [] + + for run_id, event_file in enumerate(event_files): + + trial_key = f'gamble_run{run_id + 1}' + gain_key = f'gain_run{run_id + 1}' + loss_key = f'loss_run{run_id + 1}' + + onsets.update({trial_key: []}) + durations.update({trial_key: []}) + weights_gain.update({gain_key: []}) + weights_loss.update({loss_key: []}) + + with open(event_file, 'rt') as file: + next(file) # skip the header + + for line in file: + info = line.strip().split() + onsets[trial_key].append(float(info[0])) + durations[trial_key].append(float(info[1])) + weights_gain[gain_key].append(float(info[2])) + weights_loss[loss_key].append(float(info[3])) + + # Create a Bunch per run, i.e. cond1_run1, cond2_run1, etc. + subject_info.append( + Bunch( + conditions = ['gamble'], + onsets = [onsets[trial_key]], + durations = [durations[trial_key]], + amplitudes = None, + tmod = None, + pmod = [Bunch( + name = ['gain', 'loss'], + poly = [1, 1], + param = [weights_gain[gain_key], weights_loss[loss_key]] + )], + 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 + + # 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'] + ]))) + + # 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') + + # 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(gunzip, 'out_file', specify_model, 'functional_runs') + + model_design = Node(Level1Design(), name = 'model_design') + model_design.inputs.bases = {'hrf': {'derivs': [0, 0]}} + 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') + + smooth = Node(Smooth(), name = 'smooth') + smooth.inputs.fwhm = self.fwhm + smooth.overwrite = False + subject_level.connect(contrast_estimate, 'con_images', smooth, 'in_files') + + # 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( + smooth, 'smoothed_files', 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'] + ) + + # Add connections + subject_level.connect([ + (data_sink, remove_gunzip, [('out_file', '_')]), + (gunzip, remove_gunzip, [('out_file', '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'scon_{contrast_id}.nii')\ + for contrast_id in self.contrast_list] + templates += [join( + self.directories.output_dir, + 'subject_level', '_subject_id_{subject_id}', 'SPM.mat')] + + # 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_*', 'scon_{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.extent_threshold = 5 + 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_*', 'scon_{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 + threshold.inputs.extent_threshold = 5 + + # 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', '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_3TR7.py b/tests/pipelines/test_team_3TR7.py new file mode 100644 index 00000000..1159ec8c --- /dev/null +++ b/tests/pipelines/test_team_3TR7.py @@ -0,0 +1,122 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Tests of the 'narps_open.pipelines.team_3TR7' module. + +Launch this test with PyTest + +Usage: +====== + pytest -q test_team_3TR7.py + pytest -q test_team_3TR7.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_3TR7 import PipelineTeam3TR7 + +class TestPipelinesTeam3TR7: + """ A class that contains all the unit tests for the PipelineTeam3TR7 class.""" + + @staticmethod + @mark.unit_test + def test_create(): + """ Test the creation of a PipelineTeam3TR7 object """ + + pipeline = PipelineTeam3TR7() + + # 1 - check the parameters + assert pipeline.fwhm == 8.0 + assert pipeline.team_id == '3TR7' + + # 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 PipelineTeam3TR7 object """ + pipeline = PipelineTeam3TR7() + # 1 - 1 subject outputs + pipeline.subject_list = ['001'] + helpers.test_pipeline_outputs(pipeline, [0, 0, 3, 8*2*2 + 5*2, 18]) + + # 2 - 4 subjects outputs + pipeline.subject_list = ['001', '002', '003', '004'] + helpers.test_pipeline_outputs(pipeline, [0, 0, 12, 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 = PipelineTeam3TR7.get_subject_information([test_file, test_file]) + + # Compare bunches to expected + bunch = info[0] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['gamble'] + 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', 'loss'] + assert bunch.pmod[0].poly == [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]]) + assert bunch.regressor_names is None + assert bunch.regressors is None + + bunch = info[1] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['gamble'] + 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', 'loss'] + assert bunch.pmod[0].poly == [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]]) + 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_3TR7', 'confounds.tsv') + + # Get new confounds file + PipelineTeam3TR7.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 PipelineTeam3TR7 and compare results """ + helpers.test_pipeline_evaluation('3TR7') diff --git a/tests/test_data/pipelines/team_3TR7/confounds.tsv b/tests/test_data/pipelines/team_3TR7/confounds.tsv new file mode 100644 index 00000000..cf63c178 --- /dev/null +++ b/tests/test_data/pipelines/team_3TR7/confounds.tsv @@ -0,0 +1,3 @@ +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 +-2.56954e-05 -0.00923735 0.0549667 0.000997278 -0.00019745 -0.000398988 From 1ae9ce16e5125b57583ccf0272918c918a3e78f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= <117362283+bclenet@users.noreply.github.com> Date: Fri, 12 Apr 2024 09:16:26 +0200 Subject: [PATCH 5/8] L7J7 reproduction (#189) * Starting repro + testing for L7J7 * Pipeline implemented * Extent thresholds for clustering * COntrasts naming * Hypo names * Hypo names * Er vs Ei in group comp * Output file names --- narps_open/pipelines/__init__.py | 2 +- narps_open/pipelines/team_L7J7.py | 672 ++++++++++++++++++ tests/pipelines/test_team_L7J7.py | 122 ++++ .../pipelines/team_L7J7/confounds.tsv | 3 + 4 files changed, 798 insertions(+), 1 deletion(-) create mode 100644 narps_open/pipelines/team_L7J7.py create mode 100644 tests/pipelines/test_team_L7J7.py create mode 100644 tests/test_data/pipelines/team_L7J7/confounds.tsv diff --git a/narps_open/pipelines/__init__.py b/narps_open/pipelines/__init__.py index 4294853c..87a8dc70 100644 --- a/narps_open/pipelines/__init__.py +++ b/narps_open/pipelines/__init__.py @@ -55,7 +55,7 @@ 'K9P0': None, 'L1A8': None, 'L3V8': None, - 'L7J7': None, + 'L7J7': 'PipelineTeamL7J7', 'L9G5': None, 'O03M': None, 'O21U': None, diff --git a/narps_open/pipelines/team_L7J7.py b/narps_open/pipelines/team_L7J7.py new file mode 100644 index 00000000..973384f2 --- /dev/null +++ b/narps_open/pipelines/team_L7J7.py @@ -0,0 +1,672 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Write the work of NARPS' team L7J7 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 PipelineTeamL7J7(Pipeline): + """ A class that defines the pipeline of team L7J7. """ + + def __init__(self): + super().__init__() + self.fwhm = 6.0 + self.team_id = 'L7J7' + self.contrast_list = ['0001', '0002'] + self.subject_level_contrasts = [ + ('effect_of_gain', 'T', ['gamble', 'gamblexgain^1', 'gamblexloss^1'], [0, 1, 0]), + ('effect_of_loss', 'T', ['gamble', 'gamblexgain^1', 'gamblexloss^1'], [0, 0, 1]) + ] + + def get_preprocessing(self): + """ No preprocessing has been done by team L7J7 """ + return None + + def get_run_level_analysis(self): + """ No run level analysis has been done by team L7J7 """ + 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 + + onsets = {} + durations = {} + weights_gain = {} + weights_loss = {} + + subject_info = [] + + for run_id, event_file in enumerate(event_files): + + trial_key = f'gamble_run{run_id + 1}' + gain_key = f'gain_run{run_id + 1}' + loss_key = f'loss_run{run_id + 1}' + + onsets.update({trial_key: []}) + durations.update({trial_key: []}) + weights_gain.update({gain_key: []}) + weights_loss.update({loss_key: []}) + + with open(event_file, 'rt') as file: + next(file) # skip the header + + for line in file: + info = line.strip().split() + onsets[trial_key].append(float(info[0])) + durations[trial_key].append(float(info[1])) + weights_gain[gain_key].append(float(info[2])) + weights_loss[loss_key].append(float(info[3])) + + # Create a Bunch per run, i.e. cond1_run1, cond2_run1, etc. + subject_info.append( + Bunch( + conditions = ['gamble'], + onsets = [onsets[trial_key]], + durations = [durations[trial_key]], + amplitudes = None, + tmod = None, + pmod = [Bunch( + name = ['gain', 'loss'], + poly = [1, 1], + param = [weights_gain[gain_key], weights_loss[loss_key]] + )], + 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 + + # 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'] + ]))) + + # 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': [0, 0]}} + 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')] + + # 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', '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_0001_thr.nii'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0002', 'spmT_0001.nii'), + # Hypothesis 6 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0002', '_threshold1', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0002', 'spmT_0001.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_L7J7.py b/tests/pipelines/test_team_L7J7.py new file mode 100644 index 00000000..e84e32c2 --- /dev/null +++ b/tests/pipelines/test_team_L7J7.py @@ -0,0 +1,122 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Tests of the 'narps_open.pipelines.team_L7J7' module. + +Launch this test with PyTest + +Usage: +====== + pytest -q test_team_L7J7.py + pytest -q test_team_L7J7.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_L7J7 import PipelineTeamL7J7 + +class TestPipelinesTeamL7J7: + """ A class that contains all the unit tests for the PipelineTeamL7J7 class.""" + + @staticmethod + @mark.unit_test + def test_create(): + """ Test the creation of a PipelineTeamL7J7 object """ + + pipeline = PipelineTeamL7J7() + + # 1 - check the parameters + assert pipeline.fwhm == 6.0 + assert pipeline.team_id == 'L7J7' + + # 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 PipelineTeamL7J7 object """ + pipeline = PipelineTeamL7J7() + # 1 - 1 subject outputs + pipeline.subject_list = ['001'] + helpers.test_pipeline_outputs(pipeline, [0, 0, 3, 8*2*2 + 5*2, 18]) + + # 2 - 4 subjects outputs + pipeline.subject_list = ['001', '002', '003', '004'] + helpers.test_pipeline_outputs(pipeline, [0, 0, 12, 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 = PipelineTeamL7J7.get_subject_information([test_file, test_file]) + + # Compare bunches to expected + bunch = info[0] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['gamble'] + 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', 'loss'] + assert bunch.pmod[0].poly == [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]]) + assert bunch.regressor_names is None + assert bunch.regressors is None + + bunch = info[1] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['gamble'] + 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', 'loss'] + assert bunch.pmod[0].poly == [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]]) + 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_L7J7', 'confounds.tsv') + + # Get new confounds file + PipelineTeamL7J7.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 PipelineTeamL7J7 and compare results """ + helpers.test_pipeline_evaluation('L7J7') diff --git a/tests/test_data/pipelines/team_L7J7/confounds.tsv b/tests/test_data/pipelines/team_L7J7/confounds.tsv new file mode 100644 index 00000000..cf63c178 --- /dev/null +++ b/tests/test_data/pipelines/team_L7J7/confounds.tsv @@ -0,0 +1,3 @@ +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 +-2.56954e-05 -0.00923735 0.0549667 0.000997278 -0.00019745 -0.000398988 From 695668a9d72e782b7062f67b3efaa1406e9acf0c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= <117362283+bclenet@users.noreply.github.com> Date: Tue, 16 Apr 2024 16:45:37 +0200 Subject: [PATCH 6/8] Adding participants exclusions in narps_open_runner (#194) * Adding a command line tool showing the correlation results of a pipeline execution * [DOC] install doc about correlation command line tool [skip ci] * Modifications on runner * Correlation main + exclusions in runner --- INSTALL.md | 6 +++ narps_open/runner.py | 16 +++++- .../__init__.py} | 0 narps_open/utils/correlation/__main__.py | 53 +++++++++++++++++++ setup.py | 1 + tests/conftest.py | 18 +++---- 6 files changed, 84 insertions(+), 10 deletions(-) rename narps_open/utils/{correlation.py => correlation/__init__.py} (100%) create mode 100644 narps_open/utils/correlation/__main__.py diff --git a/INSTALL.md b/INSTALL.md index e9f124ba..28936287 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -95,6 +95,7 @@ Finally, you are able to use the scripts of the project : * `narps_open_runner`: run pipelines * `narps_open_tester`: run a pipeline and test its results against original ones from the team +* `narps_open_correlations`: compute and display correlation between results and original ones from the team * `narps_description`: get the textual description made by a team * `narps_results`: download the original results from teams * `narps_open_status`: get status information about the development process of the pipelines @@ -107,6 +108,10 @@ narps_open_runner -t 2T6S -n 40 # and produces a report with correlation values. narps_open_tester -t 08MQ +# Compute the correlation values between results of 2T6S reproduction on 60 subjects with original ones +# WARNING : 2T6S must have been previously computed with a group of 60 subjects +narps_open_correlations -t 2T6S -n 60 + # Get the description of team C88N in markdown formatting narps_description -t C88N --md @@ -121,6 +126,7 @@ narps_open_status --json > For further information about these command line tools, read the corresponding documentation pages. > * `narps_open_runner` : [docs/running.md](docs/running.md) > * `narps_open_tester` : [docs/testing.md](docs/testing.md#command-line-tool) +> * `narps_open_correlations` : [docs/correlation.md](docs/correlation.md#command-line-tool) > * `narps_description` : [docs/description.md](docs/description.md) > * `narps_results` : [docs/data.md](docs/data.md#results-from-narps-teams) > * `narps_open_status` : [docs/status.md](docs/status.md) diff --git a/narps_open/runner.py b/narps_open/runner.py index bf557ba0..597d1144 100644 --- a/narps_open/runner.py +++ b/narps_open/runner.py @@ -178,8 +178,15 @@ def main(): help='run the first levels only (preprocessing + subjects + runs)') parser.add_argument('-c', '--check', action='store_true', required=False, help='check pipeline outputs (runner is not launched)') + parser.add_argument('-e', '--exclusions', action='store_true', required=False, + help='run the analyses without the excluded subjects') arguments = parser.parse_args() + # Check arguments + if arguments.exclusions and not arguments.nsubjects: + print('Argument -e/--exclusions only works with -n/--nsubjects') + return + # Initialize a PipelineRunner runner = PipelineRunner(team_id = arguments.team) runner.pipeline.directories.dataset_dir = Configuration()['directories']['dataset'] @@ -193,7 +200,14 @@ def main(): elif arguments.rsubjects is not None: runner.random_nb_subjects = int(arguments.rsubjects) else: - runner.nb_subjects = int(arguments.nsubjects) + if arguments.exclusions: + # Intersection between the requested subset and the list of not excluded subjects + runner.subjects = list( + set(get_participants_subset(int(arguments.nsubjects))) + & set(get_participants(arguments.team)) + ) + else: + runner.nb_subjects = int(arguments.nsubjects) # Check data if arguments.check: diff --git a/narps_open/utils/correlation.py b/narps_open/utils/correlation/__init__.py similarity index 100% rename from narps_open/utils/correlation.py rename to narps_open/utils/correlation/__init__.py diff --git a/narps_open/utils/correlation/__main__.py b/narps_open/utils/correlation/__main__.py new file mode 100644 index 00000000..d086499b --- /dev/null +++ b/narps_open/utils/correlation/__main__.py @@ -0,0 +1,53 @@ +#!/usr/bin/python +# coding: utf-8 + +""" A command line tool for the narps_open.utils.correlation module """ + +from os.path import join +from argparse import ArgumentParser + +from narps_open.data.results import ResultsCollection +from narps_open.utils.configuration import Configuration +from narps_open.utils.correlation import get_correlation_coefficient +from narps_open.pipelines import get_implemented_pipelines +from narps_open.runner import PipelineRunner + +def main(): + """ Entry-point for the command line tool narps_open_correlations """ + + # Parse arguments + parser = ArgumentParser(description = 'Compare reproduced files to original results.') + parser.add_argument('-t', '--team', type = str, required = True, + help = 'the team ID', choices = get_implemented_pipelines()) + subjects.add_argument('-n', '--nsubjects', type=str, required = True, + help='the number of subjects to be selected') + arguments = parser.parse_args() + + # Initialize pipeline + runner = PipelineRunner(arguments.team) + runner.pipeline.directories.dataset_dir = Configuration()['directories']['dataset'] + runner.pipeline.directories.results_dir = Configuration()['directories']['reproduced_results'] + runner.pipeline.directories.set_output_dir_with_team_id(arguments.team) + runner.pipeline.directories.set_working_dir_with_team_id(arguments.team) + runner.nb_subjects = arguments.nsubjects + + # Indices and keys to the unthresholded maps + indices = list(range(1, 18, 2)) + + # Retrieve the paths to the reproduced files + reproduced_files = runner.pipeline.get_hypotheses_outputs() + reproduced_files = [reproduced_files[i] for i in indices] + + # Retrieve the paths to the results files + collection = ResultsCollection(arguments.team) + file_keys = [f'hypo{h}_unthresh.nii.gz' for h in range(1,10)] + results_files = [join(collection.directory, k) for k in file_keys] + + # Compute the correlation coefficients + print([ + get_correlation_coefficient(reproduced_file, results_file) + for reproduced_file, results_file in zip(reproduced_files, results_files) + ]) + +if __name__ == '__main__': + main() diff --git a/setup.py b/setup.py index b17409b6..e3c65bb0 100644 --- a/setup.py +++ b/setup.py @@ -71,6 +71,7 @@ 'narps_open_runner = narps_open.runner:main', 'narps_open_tester = narps_open.tester:main', 'narps_open_status = narps_open.utils.status:main', + 'narps_open_correlations = narps_open.utils.correlation.__main__:main', 'narps_description = narps_open.data.description.__main__:main', 'narps_results = narps_open.data.results.__main__:main' ] diff --git a/tests/conftest.py b/tests/conftest.py index f12f77a0..3e5570ff 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -22,6 +22,7 @@ from narps_open.utils.correlation import get_correlation_coefficient from narps_open.utils.configuration import Configuration from narps_open.data.results import ResultsCollection +from narps_open.data.participants import get_participants_subset # Init configuration, to ensure it is in testing mode Configuration(config_type='testing') @@ -88,13 +89,12 @@ def test_pipeline_execution( TODO : how to keep intermediate files of the low level for the next numbers of subjects ? - keep intermediate levels : boolean in PipelineRunner """ - # A list of number of subject to iterate over - nb_subjects_list = list(range( - Configuration()['testing']['pipelines']['nb_subjects_per_group'], - nb_subjects, - Configuration()['testing']['pipelines']['nb_subjects_per_group']) - ) - nb_subjects_list.append(nb_subjects) + # Create subdivisions of the requested subject list + nb_subjects_per_group = Configuration()['testing']['pipelines']['nb_subjects_per_group'] + all_subjects = get_participants_subset(nb_subjects) + subjects_lists = [] + for index in range(0, len(all_subjects), nb_subjects_per_group): + subjects_lists.append(all_subjects[index:index+nb_subjects_per_group]) # Initialize the pipeline runner = PipelineRunner(team_id) @@ -104,8 +104,8 @@ def test_pipeline_execution( runner.pipeline.directories.set_working_dir_with_team_id(team_id) # Run first level by (small) sub-groups of subjects - for subjects in nb_subjects_list: - runner.nb_subjects = subjects + for subjects_list in subjects_lists: + runner.subjects = subjects_list # Run as long as there are missing files after first level (with a max number of trials) # TODO : this is a workaround From 57fea46e1e43f31302c923dc7bc15c6bad4714a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= <117362283+bclenet@users.noreply.github.com> Date: Fri, 19 Apr 2024 10:06:09 +0200 Subject: [PATCH 7/8] O6R6 reproduction (#192) * First version of O6R6 [skip ci] * Group level of O6R6 * Bug in subject level analysis * Mask input in subject level * Input mask names in subject level analysis * Changes on group level analysis for O6R6 * Remove useless smoothing * Tstat outputs of randomise [skip ci] --- narps_open/pipelines/__init__.py | 2 +- narps_open/pipelines/team_O6R6.py | 727 +++++++++++++++++++++++ narps_open/utils/correlation/__main__.py | 2 +- tests/pipelines/test_team_O6R6.py | 116 ++++ 4 files changed, 845 insertions(+), 2 deletions(-) create mode 100644 narps_open/pipelines/team_O6R6.py create mode 100644 tests/pipelines/test_team_O6R6.py diff --git a/narps_open/pipelines/__init__.py b/narps_open/pipelines/__init__.py index 87a8dc70..348337a2 100644 --- a/narps_open/pipelines/__init__.py +++ b/narps_open/pipelines/__init__.py @@ -59,7 +59,7 @@ 'L9G5': None, 'O03M': None, 'O21U': None, - 'O6R6': None, + 'O6R6': 'PipelineTeamO6R6', 'P5F3': None, 'Q58J': None, 'Q6O0': 'PipelineTeamQ6O0', diff --git a/narps_open/pipelines/team_O6R6.py b/narps_open/pipelines/team_O6R6.py new file mode 100644 index 00000000..ec78a7d8 --- /dev/null +++ b/narps_open/pipelines/team_O6R6.py @@ -0,0 +1,727 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Write the work of NARPS team O6R6 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, + 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 PipelineTeamO6R6(Pipeline): + """ A class that defines the pipeline of team O6R6 """ + + def __init__(self): + super().__init__() + self.team_id = 'O6R6' + self.contrast_list = ['1', '2'] + self.run_level_contrasts = [ + ('effect_of_gain', 'T', ['gain_trial', 'loss_trial'], [1, 0]), + ('effect_of_loss', 'T', ['gain_trial', 'loss_trial'], [0, 1]) + ] + + def get_preprocessing(self): + """ No preprocessing has been done by team O6R6 """ + return None + + def get_subject_information(event_file, group): + """ + Create Bunchs for specifyModel. + + Parameters : + - event_file : str, file corresponding to the run and the subject to analyze + - group : str, the group to which belong the subject + + Returns : + - subject_info : list of Bunch for 1st level analysis. + """ + from nipype.interfaces.base import Bunch + + onsets_trial = [] + durations_trial = [] + weights_trial = [] + onsets_gain_trial = [] + durations_gain_trial = [] + weights_gain_trial = [] + onsets_loss_trial = [] + durations_loss_trial = [] + weights_loss_trial = [] + + with open(event_file, 'rt') as file: + next(file) # skip the header + + for line in file: + info = line.strip().split() + + if float(info[4]) > 0.0: # Response time exists + onsets_trial.append(float(info[0])) + durations_trial.append(float(info[4])) + weights_trial.append(1.0) + + gain_amount = float(info[2]) + loss_amount = float(info[3]) + gain_weight = 1.0 + loss_weight = 1.0 + if group == 'equalIndifference': + gain_weight = (gain_amount/2) - 4.0 + loss_weight = (loss_amount/2) - 4.0 + elif group == 'equalRange': + gain_weight = gain_amount - 4.0 + loss_weight = loss_amount - 4.0 + + # Classify trial as "gain" or "loss" + if (gain_amount / loss_amount) > 1.93: # "Gain" trial + onsets_gain_trial.append(float(info[0])) + durations_gain_trial.append(float(info[4])) + weights_gain_trial.append(gain_weight) + else: # "Loss" trial + onsets_loss_trial.append(float(info[0])) + durations_loss_trial.append(float(info[4])) + weights_loss_trial.append(loss_weight) + + return [ + Bunch( + conditions = ['trial', 'gain_trial', 'loss_trial'], + onsets = [onsets_trial, onsets_gain_trial, onsets_loss_trial], + durations = [durations_trial, durations_gain_trial, durations_loss_trial], + amplitudes = [weights_trial, weights_gain_trial, weights_loss_trial] + ) + ] + + def get_subject_group(subject_id: str): + """ + Return the group of the subject (either 'equalRange' or 'equalIndifference'). + + Parameters : + - subject_id : str, the subject identifier + + Returns : + - group : str, the group to which belong the subject + """ + from narps_open.data.participants import get_group + + if subject_id in get_group('equalRange'): + return 'equalRange' + return 'equalIndifference' + + 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') + } + 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') + + # Function Node get_subject_group + # This returns the name of the subject's group + subject_group = Node(Function( + function = self.get_subject_group, + input_names = ['subject_id'], + output_names = ['group'] + ), + name = 'subject_group' + ) + run_level.connect(information_source, 'subject_id', subject_group, 'subject_id') + + # Get Subject Info - get subject specific condition information + subject_information = Node(Function( + function = self.get_subject_information, + input_names = ['event_file', 'group'], + output_names = ['subject_info'] + ), name = 'subject_information') + run_level.connect(select_files, 'events', subject_information, 'event_file') + run_level.connect(subject_group, 'group', subject_information, 'group') + + # 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(select_files, 'func', 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' : True }} + 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(select_files, 'func', 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') + + 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-*_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory= self.directories.dataset_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_brainmask_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_brainmask_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, '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') + + # Randomise Node - Perform clustering on statistical output + randomise = Node(Randomise(), name = 'randomise') + randomise.inputs.tfce = True + randomise.inputs.num_perm = 5000 + 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') + + # 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, + 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,'tstat_files', data_sink, + f'group_level_analysis_{method}_nsub_{nb_subjects}.@tstat_files') + + 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 + 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 + regressors_one_sample = Node( + Function( + function = self.get_one_sample_t_test_regressors, + input_names = ['subject_list'], + output_names = ['regressors'] + ), + name = 'regressors_one_sample', + ) + 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': + + # 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.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. """ + + # Handle equalRange and equalIndifference + parameters = { + 'contrast_id': self.contrast_list, + 'method': ['equalRange', 'equalIndifference'], + 'file': [ + 'randomise_tfce_corrp_tstat1.nii.gz', + 'randomise_tfce_corrp_tstat2.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': [ + 'randomise_tfce_corrp_tstat1.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', 'randomise_tfce_corrp_tstat1'), + 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', 'randomise_tfce_corrp_tstat1'), + 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', 'randomise_tfce_corrp_tstat1'), + 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', 'randomise_tfce_corrp_tstat1'), + 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', 'randomise_tfce_corrp_tstat2'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', 'zstat2.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', 'randomise_tfce_corrp_tstat2'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', 'zstat2.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', 'randomise_tfce_corrp_tstat1'), + 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', 'randomise_tfce_corrp_tstat1'), + 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', 'randomise_tfce_corrp_tstat1'), + 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/narps_open/utils/correlation/__main__.py b/narps_open/utils/correlation/__main__.py index d086499b..36ff4aa0 100644 --- a/narps_open/utils/correlation/__main__.py +++ b/narps_open/utils/correlation/__main__.py @@ -19,7 +19,7 @@ def main(): parser = ArgumentParser(description = 'Compare reproduced files to original results.') parser.add_argument('-t', '--team', type = str, required = True, help = 'the team ID', choices = get_implemented_pipelines()) - subjects.add_argument('-n', '--nsubjects', type=str, required = True, + parser.add_argument('-n', '--nsubjects', type=int, required = True, help='the number of subjects to be selected') arguments = parser.parse_args() diff --git a/tests/pipelines/test_team_O6R6.py b/tests/pipelines/test_team_O6R6.py new file mode 100644 index 00000000..ef13926b --- /dev/null +++ b/tests/pipelines/test_team_O6R6.py @@ -0,0 +1,116 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Tests of the 'narps_open.pipelines.team_O6R6' module. + +Launch this test with PyTest + +Usage: +====== + pytest -q test_team_O6R6.py + pytest -q test_team_O6R6.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_O6R6 import PipelineTeamO6R6 + +class TestPipelinesTeamO6R6: + """ A class that contains all the unit tests for the PipelineTeamO6R6 class.""" + + @staticmethod + @mark.unit_test + def test_create(): + """ Test the creation of a PipelineTeamO6R6 object """ + + pipeline = PipelineTeamO6R6() + + # 1 - check the parameters + assert pipeline.team_id == 'O6R6' + + # 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 PipelineTeamO6R6 object """ + + pipeline = PipelineTeamO6R6() + + # 1 - 1 subject outputs + pipeline.subject_list = ['001'] + helpers.test_pipeline_outputs(pipeline, [0, 2*1*4*4, 2*4*1 + 2*1, 8*4*2 + 4*4, 18]) + + # 2 - 4 subjects outputs + pipeline.subject_list = ['001', '002', '003', '004'] + helpers.test_pipeline_outputs(pipeline, [0, 2*4*4*4, 2*4*4 + 2*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') + + # Prepare several scenarii + info_ei = PipelineTeamO6R6.get_subject_information(test_file, 'equalIndifference') + info_er = PipelineTeamO6R6.get_subject_information(test_file, 'equalRange') + + # Compare bunches to expected + bunch = info_ei[0] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['trial', 'gain_trial', 'loss_trial'] + helpers.compare_float_2d_arrays(bunch.onsets, [ + [4.071, 11.834, 27.535, 36.435], + [4.071, 11.834], + [27.535, 36.435] + ]) + helpers.compare_float_2d_arrays(bunch.durations, [ + [2.388, 2.289, 2.08, 2.288], + [2.388, 2.289], + [2.08, 2.288] + ]) + helpers.compare_float_2d_arrays(bunch.amplitudes, [ + [1.0, 1.0, 1.0, 1.0], + [3.0, 13.0], + [3.5, 4.5] + ]) + + # Compare bunches to expected + bunch = info_er[0] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['trial', 'gain_trial', 'loss_trial'] + helpers.compare_float_2d_arrays(bunch.onsets, [ + [4.071, 11.834, 27.535, 36.435], + [4.071, 11.834], + [27.535, 36.435] + ]) + helpers.compare_float_2d_arrays(bunch.durations, [ + [2.388, 2.289, 2.08, 2.288], + [2.388, 2.289], + [2.08, 2.288] + ]) + helpers.compare_float_2d_arrays(bunch.amplitudes, [ + [1.0, 1.0, 1.0, 1.0], + [10.0, 30.0], + [11.0, 13.0] + ]) + + @staticmethod + @mark.pipeline_test + def test_execution(): + """ Test the execution of a PipelineTeamO6R6 and compare results """ + helpers.test_pipeline_evaluation('O6R6') From f35ea9720f7b86b2e9f3dafdf40271609c23e0fc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Boris=20Cl=C3=A9net?= <117362283+bclenet@users.noreply.github.com> Date: Mon, 29 Apr 2024 16:04:13 +0200 Subject: [PATCH 8/8] O21U reproduction (#202) * Starting O21U reproduction [skip ci] * Group level outputs * Hypotheses names * Changing unit test for O6R6 * Changing unit test for o21U --- narps_open/pipelines/__init__.py | 2 +- narps_open/pipelines/team_O21U.py | 742 ++++++++++++++++++ tests/pipelines/test_team_O21U.py | 128 +++ tests/pipelines/test_team_O6R6.py | 4 +- .../pipelines/team_O21U/confounds.tsv | 3 + 5 files changed, 876 insertions(+), 3 deletions(-) create mode 100644 narps_open/pipelines/team_O21U.py create mode 100644 tests/pipelines/test_team_O21U.py create mode 100644 tests/test_data/pipelines/team_O21U/confounds.tsv diff --git a/narps_open/pipelines/__init__.py b/narps_open/pipelines/__init__.py index 348337a2..dc1b254f 100644 --- a/narps_open/pipelines/__init__.py +++ b/narps_open/pipelines/__init__.py @@ -58,7 +58,7 @@ 'L7J7': 'PipelineTeamL7J7', 'L9G5': None, 'O03M': None, - 'O21U': None, + 'O21U': 'PipelineTeamO21U', 'O6R6': 'PipelineTeamO6R6', 'P5F3': None, 'Q58J': None, diff --git a/narps_open/pipelines/team_O21U.py b/narps_open/pipelines/team_O21U.py new file mode 100644 index 00000000..6fa70c53 --- /dev/null +++ b/narps_open/pipelines/team_O21U.py @@ -0,0 +1,742 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Write the work of NARPS team O21U 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, + FSLCommand, Cluster + ) +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 PipelineTeamO21U(Pipeline): + """ A class that defines the pipeline of team O21U """ + + def __init__(self): + super().__init__() + self.fwhm = 5.0 + self.team_id = 'O21U' + self.contrast_list = ['1', '2'] + self.run_level_contrasts = [ + ('effect_of_gain', 'T', ['gain', 'loss'], [1, 0]), + ('effect_of_loss', 'T', ['gain', 'loss'], [0, 1]) + ] + + def get_preprocessing(self): + """ No preprocessing has been done by team O21U """ + 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 = [] + durations = [] + amplitudes_trial = [] + amplitudes_gain = [] + amplitudes_loss = [] + + 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])) + amplitudes_trial.append(1.0) + amplitudes_gain.append(float(info[2])) + amplitudes_loss.append(float(info[3])) + + return [ + Bunch( + conditions = ['trial', 'gain', 'loss'], + onsets = [onsets] * 3, + durations = [durations] * 3, + amplitudes = [amplitudes_trial, amplitudes_gain, amplitudes_loss] + ) + ] + + 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 + + data_frame = read_csv(filepath, sep = '\t', header=0) + retained_confounds = DataFrame(transpose(array([ + data_frame['FramewiseDisplacement'], data_frame['X'], data_frame['Y'], data_frame['Z'], + data_frame['RotX'], data_frame['RotY'], data_frame['RotZ'] + ]))) + + # 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'), + 'mask' : join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-{run_id}_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz'), + 'events' : join('sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-{run_id}_events.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') + + # MultiImageMaths Node - Apply mask to func + mask_func = Node(MultiImageMaths(), name = 'mask_func') + mask_func.inputs.op_string = '-mul %s ' + run_level.connect(select_files, 'func', mask_func, 'in_file') + run_level.connect(select_files, 'mask', mask_func, 'operand_files') + + # IsotropicSmooth Node - Smoothing data + smoothing_func = Node(IsotropicSmooth(), name = 'smoothing_func') + smoothing_func.inputs.fwhm = self.fwhm + run_level.connect(mask_func, 'out_file', 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') + + # 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(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' : True }} + 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_masked = Node( + InterfaceFactory.create('remove_parent_directory'), + name = 'remove_masked') + run_level.connect(smoothing_func, 'out_file', remove_masked, '_') + run_level.connect(mask_func, 'out_file', remove_masked, 'file_name') + + 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-*_bold_space-MNI152NLin2009cAsym_brainmask.nii.gz') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory= self.directories.dataset_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_brainmask_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_brainmask_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, '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') + + # Cluster Node - Perform clustering on statistical output + cluster = MapNode( + Cluster(), + name = 'cluster', + iterfield = ['in_file', 'cope_file'], + synchronize = True + ) + cluster.inputs.threshold = 2.3 + cluster.inputs.out_threshold_file = True + group_level.connect(estimate_model, 'zstats', cluster, 'in_file') + group_level.connect(estimate_model, 'copes', cluster, 'cope_file') + + # 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, + 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(cluster,'threshold_file', data_sink, + f'group_level_analysis_{method}_nsub_{nb_subjects}.@threshold_file') + + 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 + 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 + regressors_one_sample = Node( + Function( + function = self.get_one_sample_t_test_regressors, + input_names = ['subject_list'], + output_names = ['regressors'] + ), + name = 'regressors_one_sample', + ) + 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': + + # 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.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. """ + + # Handle equalRange and equalIndifference + parameters = { + 'contrast_id': self.contrast_list, + 'method': ['equalRange', 'equalIndifference'], + 'file': [ + '_cluster0/zstat1_threshold.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_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 = [ + # Hypothesis 1 + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', 'randomise_tfce_corrp_tstat2'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + # Hypothesis 2 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', 'randomise_tfce_corrp_tstat2'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + # Hypothesis 3 + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', 'randomise_tfce_corrp_tstat2'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + # Hypothesis 4 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', 'randomise_tfce_corrp_tstat2'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + # Hypothesis 5 + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', 'randomise_tfce_corrp_tstat1'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', 'zstat2.nii.gz'), + # Hypothesis 6 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', 'randomise_tfce_corrp_tstat1'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', 'zstat2.nii.gz'), + # Hypothesis 7 + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', 'randomise_tfce_corrp_tstat1'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz'), + # Hypothesis 8 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', 'randomise_tfce_corrp_tstat1'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz'), + # Hypothesis 9 + join(f'group_level_analysis_groupComp_nsub_{nb_sub}', + '_contrast_id_2', 'randomise_tfce_corrp_tstat1'), + 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_O21U.py b/tests/pipelines/test_team_O21U.py new file mode 100644 index 00000000..3d65297a --- /dev/null +++ b/tests/pipelines/test_team_O21U.py @@ -0,0 +1,128 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Tests of the 'narps_open.pipelines.team_O21U' module. + +Launch this test with PyTest + +Usage: +====== + pytest -q test_team_O21U.py + pytest -q test_team_O21U.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_O21U import PipelineTeamO21U + +class TestPipelinesTeamO21U: + """ A class that contains all the unit tests for the PipelineTeamO21U class.""" + + @staticmethod + @mark.unit_test + def test_create(): + """ Test the creation of a PipelineTeamO21U object """ + + pipeline = PipelineTeamO21U() + + # 1 - check the parameters + assert pipeline.fwhm == 5.0 + assert pipeline.team_id == 'O21U' + + # 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 PipelineTeamO21U object """ + + pipeline = PipelineTeamO21U() + + # 1 - 1 subject outputs + pipeline.subject_list = ['001'] + helpers.test_pipeline_outputs(pipeline, [0, 4*1*2*4, 4*2*1 + 2*1, 6*2*2 + 3*2, 18]) + + # 2 - 4 subjects outputs + pipeline.subject_list = ['001', '002', '003', '004'] + helpers.test_pipeline_outputs(pipeline, [0, 4*4*2*4, 4*2*4 + 2*4, 6*2*2 + 3*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') + + # Prepare several scenarii + info_missed = PipelineTeamO21U.get_subject_information(test_file) + + # Compare bunches to expected + bunch = info_missed[0] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['trial', 'gain', 'loss'] + helpers.compare_float_2d_arrays(bunch.onsets, [ + [4.071, 11.834, 19.535, 27.535, 36.435], + [4.071, 11.834, 19.535, 27.535, 36.435], + [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], + [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, 1.0], + [14.0, 34.0, 38.0, 10.0, 16.0], + [6.0, 14.0, 19.0, 15.0, 17.0] + ]) + + @staticmethod + @mark.unit_test + def test_confounds_file(temporary_data_dir): + """ Test the get_confounds_file method """ + + # Get input and reference output file + confounds_file = abspath(join( + Configuration()['directories']['test_data'], 'pipelines', 'confounds.tsv')) + reference_file = abspath(join( + Configuration()['directories']['test_data'], + 'pipelines', 'team_O21U', 'confounds.tsv')) + + # Create new confounds file + confounds_node = Node(Function( + input_names = ['filepath', 'subject_id', 'run_id'], + output_names = ['confounds_file'], + function = PipelineTeamO21U.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) + + # Check contents + assert cmp(reference_file, created_confounds_file) + + @staticmethod + @mark.pipeline_test + def test_execution(): + """ Test the execution of a PipelineTeamO21U and compare results """ + helpers.test_pipeline_evaluation('O21U') diff --git a/tests/pipelines/test_team_O6R6.py b/tests/pipelines/test_team_O6R6.py index ef13926b..bd571d4d 100644 --- a/tests/pipelines/test_team_O6R6.py +++ b/tests/pipelines/test_team_O6R6.py @@ -51,11 +51,11 @@ def test_outputs(): # 1 - 1 subject outputs pipeline.subject_list = ['001'] - helpers.test_pipeline_outputs(pipeline, [0, 2*1*4*4, 2*4*1 + 2*1, 8*4*2 + 4*4, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 2*1*4*4, 2*4*1 + 2*1, 6*2*2 + 3*2, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - helpers.test_pipeline_outputs(pipeline, [0, 2*4*4*4, 2*4*4 + 2*4, 8*4*2 + 4*4, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 2*4*4*4, 2*4*4 + 2*4, 6*2*2 + 3*2, 18]) @staticmethod @mark.unit_test diff --git a/tests/test_data/pipelines/team_O21U/confounds.tsv b/tests/test_data/pipelines/team_O21U/confounds.tsv new file mode 100644 index 00000000..04e5f158 --- /dev/null +++ b/tests/test_data/pipelines/team_O21U/confounds.tsv @@ -0,0 +1,3 @@ +0.0 0.0 0.0 0.0 0.0 -0.0 0.0 +0.1352790093099999 -0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 +0.12437666391 -2.56954e-05 -0.00923735 0.0549667 0.000997278 -0.00019745 -0.000398988