diff --git a/.github/workflows/code_quality.yml b/.github/workflows/code_quality.yml index a9248671..555396c3 100644 --- a/.github/workflows/code_quality.yml +++ b/.github/workflows/code_quality.yml @@ -46,7 +46,7 @@ jobs: - name: Analyse the code with pylint run: | - pylint --fail-under 8 --ignore-paths narps_open/pipelines/ narps_open > pylint_report_narps_open.txt + pylint --fail-under 8 narps_open > pylint_report_narps_open.txt pylint --fail-under 8 tests > pylint_report_tests.txt - name: Archive pylint results diff --git a/narps_open/pipelines/__init__.py b/narps_open/pipelines/__init__.py index 936dd470..e1232851 100644 --- a/narps_open/pipelines/__init__.py +++ b/narps_open/pipelines/__init__.py @@ -69,12 +69,12 @@ 'R9K3': None, 'SM54': None, 'T54A': 'PipelineTeamT54A', - 'U26C': None, + 'U26C': 'PipelineTeamU26C', 'UI76': None, 'UK24': None, 'V55J': None, 'VG39': None, - 'X19V': None, + 'X19V': 'PipelineTeamX19V', 'X1Y5': None, 'X1Z4': None, 'XU70': None diff --git a/narps_open/pipelines/__main__.py b/narps_open/pipelines/__main__.py deleted file mode 100644 index 60fd5c76..00000000 --- a/narps_open/pipelines/__main__.py +++ /dev/null @@ -1,30 +0,0 @@ -#!/usr/bin/python -# coding: utf-8 - -""" Provide a command-line interface for the package narps_open.pipelines """ - -from argparse import ArgumentParser - -from narps_open.pipelines import get_implemented_pipelines - -def main(): - """ Entry-point for the command line tool narps_open_pipeline """ - - # Parse arguments - parser = ArgumentParser(description='Get description of a NARPS pipeline.') - parser.add_argument('-v', '--verbose', action='store_true', - help='verbose mode') - arguments = parser.parse_args() - - # Print header - print('NARPS Open Pipelines') - - # Print general information about NARS Open Pipelines - print('A codebase reproducing the 70 pipelines of the NARPS study (Botvinik-Nezer et al., 2020) shared as an open resource for the community.') - - # Print pipelines - implemented_pipelines = get_implemented_pipelines() - print(f'There are currently {len(implemented_pipelines)} implemented pipelines: {implemented_pipelines}') - -if __name__ == '__main__': - main() diff --git a/narps_open/pipelines/team_08MQ.py b/narps_open/pipelines/team_08MQ.py index 9766c3ce..0a47ca67 100644 --- a/narps_open/pipelines/team_08MQ.py +++ b/narps_open/pipelines/team_08MQ.py @@ -386,23 +386,24 @@ def get_preprocessing_outputs(self): parameters = { 'subject_id': self.subject_list, 'run_id': self.run_list, - 'file': [ - 'components_file.txt', - 'sub-{subject_id}_task-MGT_run-{run_id}_bold_brain_mcf.nii.gz.par', - 'sub-{subject_id}_task-MGT_run-{run_id}_bold_brain_mcf_st_smooth_flirt_wtsimt.nii.gz', - 'sub-{subject_id}_task-MGT_run-{run_id}_bold_brain_mask_flirt_wtsimt.nii.gz' - ] } parameter_sets = product(*parameters.values()) - template = join( + output_dir = join( self.directories.output_dir, 'preprocessing', - '_run_id_{run_id}_subject_id_{subject_id}', - '{file}' - ) + '_run_id_{run_id}_subject_id_{subject_id}' + ) + templates = [ + join(output_dir, 'components_file.txt'), + join(output_dir, 'sub-{subject_id}_task-MGT_run-{run_id}_bold_brain_mcf.nii.gz.par'), + join(output_dir, + 'sub-{subject_id}_task-MGT_run-{run_id}_bold_brain_mcf_st_smooth_flirt_wtsimt.nii.gz'), + join(output_dir, + 'sub-{subject_id}_task-MGT_run-{run_id}_bold_brain_mask_flirt_wtsimt.nii.gz') + ] return [template.format(**dict(zip(parameters.keys(), parameter_values)))\ - for parameter_values in parameter_sets] + for parameter_values in parameter_sets for template in templates] def get_subject_information(event_file): """ @@ -586,21 +587,20 @@ def get_run_level_outputs(self): 'run_id' : self.run_list, 'subject_id' : self.subject_list, 'contrast_id' : self.contrast_list, - 'file' : [ - join('results', 'cope{contrast_id}.nii.gz'), - join('results', 'tstat{contrast_id}.nii.gz'), - join('results', 'varcope{contrast_id}.nii.gz'), - join('results', 'zstat{contrast_id}.nii.gz'), - ] } parameter_sets = product(*parameters.values()) - template = join( + output_dir = join( self.directories.output_dir, - 'run_level_analysis', '_run_id_{run_id}_subject_id_{subject_id}','{file}' + '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_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ - for parameter_values in parameter_sets] + for parameter_values in parameter_sets for template in templates] return return_list @@ -1005,6 +1005,57 @@ def get_group_level_analysis_sub_workflow(self, method): 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': [ + 'randomise_tfce_corrp_tstat1.nii.gz', + 'randomise_tfce_corrp_tstat2.nii.gz', + 'randomise_tstat1.nii.gz', + 'randomise_tstat2.nii.gz', + 'tstat1.nii.gz', + 'tstat2.nii.gz', + 'zstat1.nii.gz', + 'zstat2.nii.gz' + ], + '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, + 'file': [ + 'randomise_tfce_corrp_tstat1.nii.gz', + 'randomise_tstat1.nii.gz', + 'zstat1.nii.gz', + 'tstat1.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 the names of the files used by the team to answer the hypotheses of NARPS. """ diff --git a/narps_open/pipelines/team_0I4U_debug.py b/narps_open/pipelines/team_0I4U_debug.py index bdd1b3e0..98b51837 100755 --- a/narps_open/pipelines/team_0I4U_debug.py +++ b/narps_open/pipelines/team_0I4U_debug.py @@ -1,3 +1,4 @@ +# pylint: skip-file from nipype.interfaces.spm import (Coregister, Smooth, OneSampleTTestDesign, EstimateModel, EstimateContrast, Level1Design, TwoSampleTTestDesign, RealignUnwarp, FieldMap, NewSegment, Normalize12, Reslice) diff --git a/narps_open/pipelines/team_1KB2_debug.py b/narps_open/pipelines/team_1KB2_debug.py index d56b939b..d7699667 100755 --- a/narps_open/pipelines/team_1KB2_debug.py +++ b/narps_open/pipelines/team_1KB2_debug.py @@ -1,3 +1,4 @@ +# pylint: skip-file from nipype.interfaces.fsl import (BET, FAST, MCFLIRT, FLIRT, FNIRT, ApplyWarp, SUSAN, Info, ImageMaths, IsotropicSmooth, Threshold, Level1Design, FEATModel, L2Model, Merge, FLAMEO, ContrastMgr,Cluster, FILMGLS, Randomise, MultipleRegressDesign) diff --git a/narps_open/pipelines/team_43FJ_debug.py b/narps_open/pipelines/team_43FJ_debug.py index 08be5ec9..011e4e63 100755 --- a/narps_open/pipelines/team_43FJ_debug.py +++ b/narps_open/pipelines/team_43FJ_debug.py @@ -1,3 +1,4 @@ +# pylint: skip-file from nipype.interfaces.fsl import (BET, FAST, MCFLIRT, FLIRT, FNIRT, ApplyWarp, SUSAN, MotionOutliers, Info, ImageMaths, IsotropicSmooth, Threshold, Level1Design, FEATModel, L2Model, Merge, FLAMEO, ContrastMgr, FILMGLS, Randomise, MultipleRegressDesign) diff --git a/narps_open/pipelines/team_4TQ6_wip.py b/narps_open/pipelines/team_4TQ6_wip.py index 623df025..28e9aada 100755 --- a/narps_open/pipelines/team_4TQ6_wip.py +++ b/narps_open/pipelines/team_4TQ6_wip.py @@ -1,3 +1,4 @@ +# pylint: skip-file from nipype.interfaces.fsl import (BET, ICA_AROMA, FAST, MCFLIRT, FLIRT, FNIRT, ApplyWarp, SUSAN, Info, ImageMaths, IsotropicSmooth, Threshold, Level1Design, FEATModel, L2Model, Merge, FLAMEO, ContrastMgr,Cluster, FILMGLS, Randomise, MultipleRegressDesign) diff --git a/narps_open/pipelines/team_Q6O0.py b/narps_open/pipelines/team_Q6O0.py index 69cacc3c..176b84cf 100755 --- a/narps_open/pipelines/team_Q6O0.py +++ b/narps_open/pipelines/team_Q6O0.py @@ -468,11 +468,14 @@ def get_subject_level_outputs(self): ) # Formatting templates and returning it as a list of files - output_files = [contrast_map_template.format(**dict(zip(parameters.keys(), parameter_values)))\ + output_files = [ + contrast_map_template.format(**dict(zip(parameters.keys(), parameter_values)))\ for parameter_values in parameter_sets] - output_files += [mat_file_template.format(**dict(zip(parameters.keys(), parameter_values)))\ + output_files += [ + mat_file_template.format(**dict(zip(parameters.keys(), parameter_values)))\ for parameter_values in parameter_sets] - output_files += [spmt_file_template.format(**dict(zip(parameters.keys(), parameter_values)))\ + output_files += [ + spmt_file_template.format(**dict(zip(parameters.keys(), parameter_values)))\ for parameter_values in parameter_sets] return output_files @@ -698,23 +701,41 @@ def get_hypotheses_outputs(self): """ Return all hypotheses output file names. """ nb_sub = len(self.subject_list) files = [ - join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', '_model_type_gain', '_threshold0', 'spmT_0001_thr.nii'), - join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', '_model_type_gain', 'spmT_0001.nii'), - join(f'l2_analysis_equalRange_nsub_{nb_sub}', '_model_type_gain', '_threshold0', 'spmT_0001_thr.nii'), - join(f'l2_analysis_equalRange_nsub_{nb_sub}', '_model_type_gain', 'spmT_0001.nii'), - join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', '_model_type_gain', '_threshold0', 'spmT_0001_thr.nii'), - join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', '_model_type_gain', 'spmT_0001.nii'), - join(f'l2_analysis_equalRange_nsub_{nb_sub}', '_model_type_gain', '_threshold0', 'spmT_0001_thr.nii'), - join(f'l2_analysis_equalRange_nsub_{nb_sub}', '_model_type_gain', 'spmT_0001.nii'), - join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', '_model_type_loss', '_threshold1', 'spmT_0002_thr.nii'), - join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', '_model_type_loss', 'spmT_0002.nii'), - join(f'l2_analysis_equalRange_nsub_{nb_sub}', '_model_type_loss', '_threshold1', 'spmT_0002_thr.nii'), - join(f'l2_analysis_equalRange_nsub_{nb_sub}', '_model_type_loss', 'spmT_0002.nii'), - join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', '_model_type_loss', '_threshold0', 'spmT_0001_thr.nii'), - join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', '_model_type_loss', 'spmT_0001.nii'), - join(f'l2_analysis_equalRange_nsub_{nb_sub}', '_model_type_loss', '_threshold0', 'spmT_0001_thr.nii'), - join(f'l2_analysis_equalRange_nsub_{nb_sub}', '_model_type_loss', 'spmT_0001.nii'), - join(f'l2_analysis_groupComp_nsub_{nb_sub}', '_model_type_loss', '_threshold0', 'spmT_0001_thr.nii'), - join(f'l2_analysis_groupComp_nsub_{nb_sub}', '_model_type_loss', 'spmT_0001.nii') + join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', + '_model_type_gain', '_threshold0', 'spmT_0001_thr.nii'), + join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', + '_model_type_gain', 'spmT_0001.nii'), + join(f'l2_analysis_equalRange_nsub_{nb_sub}', + '_model_type_gain', '_threshold0', 'spmT_0001_thr.nii'), + join(f'l2_analysis_equalRange_nsub_{nb_sub}', + '_model_type_gain', 'spmT_0001.nii'), + join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', + '_model_type_gain', '_threshold0', 'spmT_0001_thr.nii'), + join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', + '_model_type_gain', 'spmT_0001.nii'), + join(f'l2_analysis_equalRange_nsub_{nb_sub}', + '_model_type_gain', '_threshold0', 'spmT_0001_thr.nii'), + join(f'l2_analysis_equalRange_nsub_{nb_sub}', + '_model_type_gain', 'spmT_0001.nii'), + join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', + '_model_type_loss', '_threshold1', 'spmT_0002_thr.nii'), + join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', + '_model_type_loss', 'spmT_0002.nii'), + join(f'l2_analysis_equalRange_nsub_{nb_sub}', + '_model_type_loss', '_threshold1', 'spmT_0002_thr.nii'), + join(f'l2_analysis_equalRange_nsub_{nb_sub}', + '_model_type_loss', 'spmT_0002.nii'), + join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', + '_model_type_loss', '_threshold0', 'spmT_0001_thr.nii'), + join(f'l2_analysis_equalIndifference_nsub_{nb_sub}', + '_model_type_loss', 'spmT_0001.nii'), + join(f'l2_analysis_equalRange_nsub_{nb_sub}', + '_model_type_loss', '_threshold0', 'spmT_0001_thr.nii'), + join(f'l2_analysis_equalRange_nsub_{nb_sub}', + '_model_type_loss', 'spmT_0001.nii'), + join(f'l2_analysis_groupComp_nsub_{nb_sub}', + '_model_type_loss', '_threshold0', 'spmT_0001_thr.nii'), + join(f'l2_analysis_groupComp_nsub_{nb_sub}', + '_model_type_loss', '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 index 8294ec2b..65613e70 100755 --- a/narps_open/pipelines/team_R9K3_wip.py +++ b/narps_open/pipelines/team_R9K3_wip.py @@ -1,3 +1,4 @@ +# pylint: skip-file # THIS IS A TEMPLATE THAT CAN BE USE TO REPRODUCE A NEW PIPELINE import os diff --git a/narps_open/pipelines/team_T54A.py b/narps_open/pipelines/team_T54A.py index ae34b848..526162ca 100644 --- a/narps_open/pipelines/team_T54A.py +++ b/narps_open/pipelines/team_T54A.py @@ -765,17 +765,22 @@ def get_group_level_outputs(self): for parameter_values in parameter_sets] # Handle groupComp - files = [ - 'randomise_tfce_corrp_tstat1.nii.gz', - 'randomise_tstat1.nii.gz', - 'zstat1.nii.gz', - 'tstat1.nii.gz' + parameters = { + 'contrast_id': self.contrast_list, + 'file': [ + 'randomise_tfce_corrp_tstat1.nii.gz', + 'randomise_tstat1.nii.gz', + 'zstat1.nii.gz', + 'tstat1.nii.gz' ] - - return_list += [join( + } + parameter_sets = product(*parameters.values()) + template = join( self.directories.output_dir, f'group_level_analysis_groupComp_nsub_{len(self.subject_list)}', - '_contrast_id_2', f'{file}') for file in files] + '_contrast_id_{contrast_id}', '{file}') + return_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] return return_list diff --git a/narps_open/pipelines/team_U26C.py b/narps_open/pipelines/team_U26C.py new file mode 100755 index 00000000..4c4e8e73 --- /dev/null +++ b/narps_open/pipelines/team_U26C.py @@ -0,0 +1,692 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Write the work of NARPS' team U26C 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 PipelineTeamU26C(Pipeline): + """ A class that defines the pipeline of team U26C. """ + + def __init__(self): + super().__init__() + self.fwhm = 5.0 + self.team_id = 'U26C' + self.contrast_list = ['0001', '0002', '0003'] + + gamble = [f'gamble_run{r}' for r in range(1, len(self.run_list) + 1)] + gain = [f'gamble_run{r}xgain_run{r}^1' for r in range(1, len(self.run_list) + 1)] + loss = [f'gamble_run{r}xloss_run{r}^1' for r in range(1, len(self.run_list) + 1)] + + self.subject_level_contrasts = [ + ['gamble', 'T', gamble, [1, 1, 1, 1]], + ['gain', 'T', gain, [1, 1, 1, 1]], + ['loss', 'T', loss, [1, 1, 1, 1]] + ] + + def get_preprocessing(self): + """ No preprocessing has been done by team U26C """ + return None + + def get_run_level_analysis(self): + """ No run level analysis has been done by team U26C """ + 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_information : 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 = [trial_key], + onsets = [onsets[trial_key]], + durations = [durations[trial_key]], + amplitudes = None, + tmod = None, + pmod = [Bunch( + name = [gain_key, loss_key], + 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, insert, diff + + # 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['CSF'], data_frame['WhiteMatter'], + data_frame['X'], data_frame['Y'], data_frame['Z'], + data_frame['RotX'], data_frame['RotY'], data_frame['RotZ'], + insert(diff(data_frame['X']), 0, 0), + insert(diff(data_frame['Y']), 0, 0), + insert(diff(data_frame['Z']), 0, 0), + insert(diff(data_frame['RotX']), 0, 0), + insert(diff(data_frame['RotY']), 0, 0), + insert(diff(data_frame['RotZ']), 0, 0) + ] + ))) + + # 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_analysis : nipype.WorkFlow + """ + # Identity interface Node - to iterate over subject_id and run + infosource = Node( + IdentityInterface(fields = ['subject_id']), + name = 'infosource') + infosource.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') + } + selectderivs = Node(SelectFiles(templates), name = 'selectderivs') + selectderivs.inputs.base_directory = self.directories.dataset_dir + selectderivs.inputs.sort_filelist = True + + # 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 + + # Gunzip - gunzip files because SPM do not use .nii.gz files + gunzip = MapNode(Gunzip(), name = 'gunzip', iterfield=['in_file']) + + # Smooth warped functionals. + smooth = Node(Smooth(), name = 'smooth') + smooth.inputs.fwhm = self.fwhm + smooth.overwrite = False + + # Function node get_subject_information - get subject specific condition information + getsubinforuns = Node(Function( + function = self.get_subject_information, + input_names = ['event_files'], + output_names = ['subject_info'] + ), + name = 'getsubinforuns') + + # Function node get_confounds_file - get 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 + + modelspec = Node(SpecifySPMModel(), name = 'modelspec') + modelspec.inputs.concatenate_runs = False + modelspec.inputs.input_units = 'secs' + modelspec.inputs.output_units = 'secs' + modelspec.inputs.time_repetition = TaskInformation()['RepetitionTime'] + modelspec.inputs.high_pass_filter_cutoff = 128 + modelspec.overwrite = False + + level1design = Node(Level1Design(), name = 'level1design') + level1design.inputs.bases = {'hrf': {'derivs': [0, 0]}} + level1design.inputs.timing_units = 'secs' + level1design.inputs.interscan_interval = TaskInformation()['RepetitionTime'] + level1design.overwrite = False + + level1estimate = Node(EstimateModel(), name = 'level1estimate') + level1estimate.inputs.estimation_method = {'Classical': 1} + level1estimate.overwrite = False + + 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_analysis = Workflow( + base_dir = self.directories.working_dir, name = 'subject_level_analysis' + ) + subject_level_analysis.connect([ + (infosource, selectderivs, [('subject_id', 'subject_id')]), + (infosource, confounds, [('subject_id', 'subject_id')]), + (selectderivs, getsubinforuns, [('events', 'event_files')]), + (selectderivs, gunzip, [('func', 'in_file')]), + (selectderivs, confounds, [('confounds', 'filepath')]), + (gunzip, smooth, [('out_file', 'in_files')]), + (getsubinforuns, modelspec, [('subject_info', 'subject_info')]), + (confounds, modelspec, [('confounds_file', 'realignment_parameters')]), + (smooth, modelspec, [('smoothed_files', 'functional_runs')]), + (modelspec, level1design, [('session_info', 'session_info')]), + (level1design, level1estimate, [('spm_mat_file', 'spm_mat_file')]), + (level1estimate, contrast_estimate,[ + ('spm_mat_file', 'spm_mat_file'), + ('beta_images', 'beta_images'), + ('residual_image', 'residual_image')]), + (contrast_estimate, data_sink, [ + ('con_images', f'{subject_level_analysis.name}.@con_images'), + ('spmT_images', f'{subject_level_analysis.name}.@spmT_images'), + ('spm_mat_file', f'{subject_level_analysis.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_analysis.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_analysis + + 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_analysis', '_subject_id_{subject_id}', f'con_{contrast_id}.nii')\ + for contrast_id in self.contrast_list] + templates += [join( + self.directories.output_dir, + 'subject_level_analysis', '_subject_id_{subject_id}', 'SPM.mat')] + 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 + """ + + 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 + infosource = Node(IdentityInterface(fields=['contrast_id']), + name = 'infosource') + infosource.iterables = [('contrast_id', self.contrast_list)] + + # Select files from subject level analysis + templates = { + 'contrasts': join(self.directories.output_dir, + 'subject_level_analysis', '_subject_id_*', 'con_{contrast_id}.nii'), + #'mask': join('derivatives/fmriprep/gr_mask_tmax.nii') + } + selectderivs = Node(SelectFiles(templates), name = 'selectderivs') + selectderivs.inputs.sort_filelist = True + selectderivs.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([ + (infosource, selectderivs, [('contrast_id', 'contrast_id')]), + (selectderivs, 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') + ]), + #(selectderivs, 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 + infosource = Node(IdentityInterface(fields=['contrast_id']), + name = 'infosource') + infosource.iterables = [('contrast_id', self.contrast_list)] + + # Select files from subject level analysis + templates = { + 'contrasts': join(self.directories.output_dir, + 'subject_level_analysis', '_subject_id_*', 'con_{contrast_id}.nii'), + #'mask': join('derivatives/fmriprep/gr_mask_tmax.nii') + } + selectderivs = Node(SelectFiles(templates), name = 'selectderivs') + selectderivs.inputs.sort_filelist = True + selectderivs.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([ + (infosource, selectderivs, [('contrast_id', 'contrast_id')]), + (selectderivs, get_equal_range_contrasts, [('contrasts', 'input_str')]), + (selectderivs, 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') + ]), + #(selectderivs, twosampttest, [('mask', 'explicit_mask_file')]), + (twosampttest, level2estimate, [('spm_mat_file', 'spm_mat_file')]), + (level2estimate, level2conestimate, [ + ('spm_mat_file', 'spm_mat_file'), + ('beta_images', 'beta_images'), + ('residual_image', 'residual_image') + ]), + (level2conestimate, threshold, [ + ('spm_mat_file', 'spm_mat_file'), + ('spmT_images', 'stat_image') + ]), + (level2estimate, data_sink, [ + ('mask_image', f'{group_level_analysis.name}.@mask')]), + (level2conestimate, data_sink, [ + ('spm_mat_file', f'{group_level_analysis.name}.@spm_mat'), + ('spmT_images', f'{group_level_analysis.name}.@T'), + ('con_images', f'{group_level_analysis.name}.@con')]), + (threshold, data_sink, [ + ('thresholded_map', f'{group_level_analysis.name}.@thresh')]) + ]) + + return group_level_analysis + + def get_group_level_outputs(self): + """ Return all names for the files the group level analysis is supposed to generate. """ + + # Handle equalRange and equalIndifference + parameters = { + 'contrast_id': self.contrast_list, + 'method': ['equalRange', 'equalIndifference'], + 'file': [ + 'con_0001.nii', 'con_0002.nii', 'mask.nii', 'SPM.mat', + 'spmT_0001.nii', 'spmT_0002.nii', + join('_threshold0', 'spmT_0001_thr.nii'), join('_threshold1', 'spmT_0002_thr.nii') + ], + 'nb_subjects' : [str(len(self.subject_list))] + } + + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'group_level_analysis_{method}_nsub_{nb_subjects}', + '_contrast_id_{contrast_id}', + '{file}' + ) + return_list = [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + # Handle groupComp + parameters = { + 'contrast_id': self.contrast_list, + 'method': ['groupComp'], + 'file': [ + 'con_0001.nii', 'mask.nii', 'SPM.mat', 'spmT_0001.nii', + join('_threshold0', 'spmT_0001_thr.nii') + ], + 'nb_subjects' : [str(len(self.subject_list))] + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'group_level_analysis_{method}_nsub_{nb_subjects}', + '_contrast_id_{contrast_id}', + '{file}' + ) + return_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + return return_list + + def get_hypotheses_outputs(self): + """ Return all hypotheses output file names. """ + nb_sub = len(self.subject_list) + files = [ + # Hypothesis 1 + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0002', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_0002', 'spmT_0001.nii'), + # Hypothesis 2 + 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 3 + 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 4 + 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 5 + 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'), + # Hypothesis 6 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0003', '_threshold1', 'spmT_0002_thr.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0003', 'spmT_0002.nii'), + # Hypothesis 7 + 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'), + # Hypothesis 8 + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0003', '_threshold0', 'spmT_0001_thr.nii'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_0003', 'spmT_0001.nii'), + # Hypothesis 9 + 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/narps_open/pipelines/team_V55J.py b/narps_open/pipelines/team_V55J.py index 9f930608..c451c77d 100755 --- a/narps_open/pipelines/team_V55J.py +++ b/narps_open/pipelines/team_V55J.py @@ -1,3 +1,4 @@ +# pylint: skip-file from nipype.interfaces.spm import (Coregister, Smooth, OneSampleTTestDesign, EstimateModel, EstimateContrast, Level1Design, TwoSampleTTestDesign, RealignUnwarp, Normalize12, NewSegment, FieldMap) diff --git a/narps_open/pipelines/team_X19V.py b/narps_open/pipelines/team_X19V.py index 3e0108ef..5e3e8989 100755 --- a/narps_open/pipelines/team_X19V.py +++ b/narps_open/pipelines/team_X19V.py @@ -1,611 +1,821 @@ -from nipype.interfaces.fsl import (Info, ImageMaths, IsotropicSmooth, Threshold, Level1Design, FEATModel, - L2Model, Merge, FLAMEO, ContrastMgr, FILMGLS, MultipleRegressDesign, Cluster, BET, SmoothEstimate) -from nipype.algorithms.modelgen import SpecifyModel -from nipype.interfaces.utility import IdentityInterface, Function -from nipype.interfaces.io import SelectFiles, DataSink -from nipype.algorithms.misc import Gunzip +#!/usr/bin/python +# coding: utf-8 + +""" Write the work of NARPS team X19V using Nipype """ + +from os.path import join +from itertools import product + from nipype import Workflow, Node, MapNode -from nipype.interfaces.base import Bunch - -from os.path import join as opj -import os - -# Event-Related design, 4 second events, with parametric modulation based on amount gained or lost on each trial. -# One 'gain' regressor, and one 'loss' regressor. -# No RT modelling - -def get_session_infos(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 os.path import join as opj - from nipype.interfaces.base import Bunch - import numpy as np - - cond_names = ['trial', 'gain', 'loss'] - - onset = {} - duration = {} - amplitude = {} - - for c in cond_names: # For each condition. - onset.update({c : []}) # creates dictionary items with empty lists - duration.update({c : []}) - amplitude.update({c : []}) - - with open(event_file, 'rt') as f: - next(f) # skip the header - - for line in f: - info = line.strip().split() - # Creates list with onsets, duration and loss/gain for amplitude (FSL) - for c in cond_names: - onset[c].append(float(info[0])) - duration[c].append(float(4)) - if c == 'gain': - amplitude[c].append(float(info[2])) - elif c == 'loss': - amplitude[c].append(float(info[3])) - elif c == 'trial': - amplitude[c].append(float(1)) - - amplitude['gain'] = amplitude['gain'] - np.mean(amplitude['gain']) - amplitude['loss'] = amplitude['loss'] - np.mean(amplitude['loss']) - - subject_info = [] - - subject_info.append(Bunch(conditions=cond_names, - onsets=[onset[k] for k in cond_names], - durations=[duration[k] for k in cond_names], - amplitudes=[amplitude[k] for k in cond_names], - regressor_names=None, - regressors=None)) - - return subject_info - -def get_parameters_file(file, subject_id, run_id, result_dir, working_dir): - ''' - Create new tsv files with only desired parameters per subject per run. - - Parameters : - - filepaths : paths to subject parameters 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 - - result_dir: str, directory where results will be stored - - working_dir: str, name of the sub-directory for intermediate results - - Return : - - parameters_file : paths to new files containing only desired parameters. - ''' - import pandas as pd - import numpy as np - from os.path import join as opj - import os - - parameters_file = [] - - df = pd.read_csv(file, sep = '\t', header=0) - temp_list = np.array([df['X'], df['Y'], df['Z'], - df['RotX'], df['RotY'], df['RotZ']]) # Parameters we want to use for the model - retained_parameters = pd.DataFrame(np.transpose(temp_list)) - new_path =opj(result_dir, working_dir, 'parameters_file', f"parameters_file_sub-{subject_id}_run{run_id}.tsv") - if not os.path.isdir(opj(result_dir, working_dir, 'parameters_file')): - os.mkdir(opj(result_dir, working_dir, 'parameters_file')) - writer = open(new_path, "w") - writer.write(retained_parameters.to_csv(sep = '\t', index = False, header = False, na_rep = '0.0')) - writer.close() - - parameters_file.append(new_path) - - return parameters_file - -# Linear contrast effects: 'Gain' vs. baseline, 'Loss' vs. baseline. -def get_contrasts(subject_id): - ''' - 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: str, ID of the subject - - Returns: - - contrasts: list of tuples, list of contrasts to analyze - ''' - # list of condition names - conditions = ['trial', 'gain', 'loss'] - - # create contrasts - gain = ('gain', 'T', conditions, [0, 1, 0]) - - loss = ('loss', 'T', conditions, [0, 0, 1]) - - gain_sup = ('gain_sup_loss', 'T', conditions, [0, 1, -1]) - - loss_sup = ('loss_sup_gain', 'T', conditions, [0, -1, 1]) - - - # contrast list - contrasts = [gain, loss, gain_sup, loss_sup] - - return contrasts - - -def rm_smoothed_files(files, subject_id, run_id, result_dir, working_dir): - import shutil - from os.path import join as opj - - smooth_dir = opj(result_dir, working_dir, 'l1_analysis', f"_run_id_{run_id}_subject_id_{subject_id}", 'smooth') - - try: - shutil.rmtree(smooth_dir) - except OSError as e: - print(e) - else: - print("The directory is deleted successfully") - - - -def get_l1_analysis(subject_list, run_list, TR, fwhm, exp_dir, output_dir, working_dir, result_dir): - """ - Returns the first level analysis workflow. - - Parameters: - - exp_dir: str, directory where raw data are stored - - result_dir: str, directory where results will be stored - - working_dir: str, name of the sub-directory for intermediate results - - output_dir: str, name of the sub-directory for final results - - subject_list: list of str, list of subject for which you want to do the analysis - - run_list: list of str, list of runs for which you want to do the analysis - - fwhm: float, fwhm for smoothing step - - TR: float, time repetition used during acquisition - - Returns: - - l1_analysis : Nipype WorkFlow - """ - # Infosource Node - To iterate on subject and runs - infosource = Node(IdentityInterface(fields = ['subject_id', 'run_id']), name = 'infosource') - infosource.iterables = [('subject_id', subject_list), - ('run_id', run_list)] - - # Templates to select files node - func_file = opj('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', - 'sub-{subject_id}_task-MGT_run-{run_id}_bold_space-MNI152NLin2009cAsym_preproc.nii.gz') - - event_file = opj('sub-{subject_id}', 'func', 'sub-{subject_id}_task-MGT_run-{run_id}_events.tsv') - - param_file = opj('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', - 'sub-{subject_id}_task-MGT_run-{run_id}_bold_confounds.tsv') - - template = {'func' : func_file, 'event' : event_file, 'param' : param_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') - - ## Skullstripping - skullstrip = Node(BET(frac = 0.1, functional = True, mask = True), name = 'skullstrip') - - ## Smoothing - smooth = Node(IsotropicSmooth(fwhm = 5), name = 'smooth') - - # Node contrasts to get contrasts - contrasts = Node(Function(function=get_contrasts, - input_names=['subject_id'], - output_names=['contrasts']), - name='contrasts') - - # Get Subject Info - get subject specific condition information - get_subject_infos = Node(Function(input_names=['event_file'], - output_names=['subject_info'], - function=get_session_infos), - name='get_subject_infos') - - specify_model = Node(SpecifyModel(high_pass_filter_cutoff = 60, - input_units = 'secs', - time_repetition = TR), name = 'specify_model') - - parameters = Node(Function(function=get_parameters_file, - input_names=['file', 'subject_id', 'run_id', 'result_dir', 'working_dir'], - output_names=['parameters_file']), - name='parameters') - - parameters.inputs.result_dir = result_dir - parameters.inputs.working_dir = working_dir - - # First temporal derivatives of the two regressors were also used, - # along with temporal filtering (60 s) of all the independent variable time-series. - # No motion parameter regressors used. - - l1_design = Node(Level1Design(bases = {'dgamma':{'derivs' : False}}, - interscan_interval = TR, - model_serial_correlations = True), name = 'l1_design') - - model_generation = Node(FEATModel(), name = 'model_generation') - - model_estimate = Node(FILMGLS(), name='model_estimate') - - remove_smooth = Node(Function(input_names = ['subject_id', 'run_id', 'files', 'result_dir', 'working_dir'], - function = rm_smoothed_files), name = 'remove_smooth') - - remove_smooth.inputs.result_dir = result_dir - remove_smooth.inputs.working_dir = working_dir - - # 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'), - ('run_id', 'run_id')]), - (selectfiles, get_subject_infos, [('event', 'event_file')]), - (selectfiles, parameters, [('param', 'file')]), - (infosource, contrasts, [('subject_id', 'subject_id')]), - (infosource, parameters, [('subject_id', 'subject_id'), - ('run_id', 'run_id')]), - (selectfiles, skullstrip, [('func', 'in_file')]), - (skullstrip, smooth, [('out_file', 'in_file')]), - (parameters, specify_model, [('parameters_file', 'realignment_parameters')]), - (smooth, specify_model, [('out_file', 'functional_runs')]), - (get_subject_infos, specify_model, [('subject_info', 'subject_info')]), - (contrasts, l1_design, [('contrasts', 'contrasts')]), - (specify_model, l1_design, [('session_info', 'session_info')]), - (l1_design, model_generation, [('ev_files', 'ev_files'), ('fsf_files', 'fsf_file')]), - (smooth, model_estimate, [('out_file', 'in_file')]), - (model_generation, model_estimate, [('con_file', 'tcon_file'), - ('design_file', 'design_file')]), - (infosource, remove_smooth, [('subject_id', 'subject_id'), ('run_id', 'run_id')]), - (model_estimate, remove_smooth, [('results_dir', 'files')]), - (model_estimate, datasink, [('results_dir', 'l1_analysis.@results')]), - (model_generation, datasink, [('design_file', 'l1_analysis.@design_file'), - ('design_image', 'l1_analysis.@design_img')]), - (skullstrip, datasink, [('mask_file', 'l1_analysis.@skullstriped')]) - ]) - - return l1_analysis - -def get_l2_analysis(subject_list, contrast_list, run_list, exp_dir, output_dir, working_dir, result_dir): - """ - Returns the 2nd level of analysis workflow. - - Parameters: - - exp_dir: str, directory where raw data are stored - - result_dir: str, directory where results will be stored - - working_dir: str, name of the sub-directory for intermediate results - - output_dir: str, name of the sub-directory for final results - - subject_list: list of str, list of subject for which you want to do the preprocessing - - contrast_list: list of str, list of contrasts to analyze - - run_list: list of str, list of runs for which you want to do the analysis - - - Returns: - - l2_analysis: Nipype WorkFlow - """ - # Infosource Node - To iterate on subject and runs - infosource_2ndlevel = Node(IdentityInterface(fields=['subject_id', 'contrast_id']), name='infosource_2ndlevel') - infosource_2ndlevel.iterables = [('subject_id', subject_list), ('contrast_id', contrast_list)] - - # Templates to select files node - copes_file = opj(output_dir, 'l1_analysis', '_run_id_*_subject_id_{subject_id}', 'results', - 'cope{contrast_id}.nii.gz') - - varcopes_file = opj(output_dir, 'l1_analysis', '_run_id_*_subject_id_{subject_id}', 'results', - 'varcope{contrast_id}.nii.gz') - - mask_file = opj(output_dir, 'l1_analysis', '_run_id_*_subject_id_{subject_id}', 'sub-{subject_id}_task-MGT_run-01_bold_space-MNI152NLin2009cAsym_preproc_brain_mask.nii.gz') - - template = {'cope' : copes_file, 'varcope' : varcopes_file, 'mask':mask_file} - - # SelectFiles node - to select necessary files - selectfiles_2ndlevel = Node(SelectFiles(template, base_directory=result_dir), name = 'selectfiles_2ndlevel') - - datasink_2ndlevel = Node(DataSink(base_directory=result_dir, container=output_dir), name='datasink_2ndlevel') - - # Generate design matrix - specify_model_2ndlevel = Node(L2Model(num_copes = len(run_list)), name='l2model_2ndlevel') - - # Merge copes and varcopes files for each subject - merge_copes_2ndlevel = Node(Merge(dimension='t'), name='merge_copes_2ndlevel') - - merge_varcopes_2ndlevel = Node(Merge(dimension='t'), name='merge_varcopes_2ndlevel') - - # Second level (single-subject, mean of all four scans) analyses: Fixed effects analysis. - flame = Node(FLAMEO(run_mode = 'flame1'), - name='flameo') - - l2_analysis = Workflow(base_dir = opj(result_dir, working_dir), name = "l2_analysis") - - l2_analysis.connect([(infosource_2ndlevel, selectfiles_2ndlevel, [('subject_id', 'subject_id'), - ('contrast_id', 'contrast_id')]), - (selectfiles_2ndlevel, merge_copes_2ndlevel, [('cope', 'in_files')]), - (selectfiles_2ndlevel, merge_varcopes_2ndlevel, [('varcope', 'in_files')]), - (selectfiles_2ndlevel, flame, [('mask', 'mask_file')]), - (merge_copes_2ndlevel, flame, [('merged_file', 'cope_file')]), - (merge_varcopes_2ndlevel, flame, [('merged_file', 'var_cope_file')]), - (specify_model_2ndlevel, flame, [('design_mat', 'design_file'), - ('design_con', 't_con_file'), - ('design_grp', 'cov_split_file')]), - (flame, datasink_2ndlevel, [('zstats', 'l2_analysis.@stats'), - ('tstats', 'l2_analysis.@tstats'), - ('copes', 'l2_analysis.@copes'), - ('var_copes', 'l2_analysis.@varcopes')])]) - - return l2_analysis - -def get_subgroups_contrasts(copes, varcopes, subject_ids, participants_file): - ''' - Parameters : - - copes: original file list selected by selectfiles node - - varcopes: original file list selected by selectfiles node - - subject_ids: list of subject IDs that are analyzed - - participants_file: str, file containing participants characteristics - - This function return the file list containing only the files belonging to subject in the wanted group. - ''' - - from os.path import join as opj - - equalRange_id = [] - equalIndifference_id = [] - - subject_list = ['sub-' + str(i) for i in subject_ids] - - with open(participants_file, 'rt') as f: - next(f) # skip the header - - for line in f: +from nipype.interfaces.utility import IdentityInterface, Function, Split +from nipype.interfaces.io import SelectFiles, DataSink +from nipype.interfaces.fsl import ( + IsotropicSmooth, Level1Design, FEATModel, + L2Model, Merge, FLAMEO, FILMGLS, MultipleRegressDesign, + Cluster, BET, SmoothEstimate, FSLCommand + ) +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 PipelineTeamX19V(Pipeline): + """ A class that defines the pipeline of team X19V """ + + def __init__(self): + super().__init__() + self.fwhm = 5.0 + self.team_id = 'X19V' + self.contrast_list = ['1', '2', '3', '4'] + self.run_level_contrasts = [ + ('gain', 'T', ['trial', 'gain', 'loss'], [0, 1, 0]), + ('loss', 'T', ['trial', 'gain', 'loss'], [0, 0, 1]), + ('gain_sup_loss', 'T', ['trial', 'gain', 'loss'], [0, 1, -1]), + ('loss_sup_gain', 'T', ['trial', 'gain', 'loss'], [0, -1, 1]) + ] + + def get_preprocessing(self): + """ No preprocessing has been done by team X19V """ + 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 + + condition_names = ['trial', 'gain', 'loss'] + onsets = {} + durations = {} + amplitudes = {} + + for condition in condition_names: + # Create dictionary items with empty lists + onsets.update({condition : []}) + durations.update({condition : []}) + amplitudes.update({condition : []}) + + with open(event_file, 'rt') as file: + next(file) # skip the header + + for line in file: info = line.strip().split() - - if info[0] in subject_list and info[1] == "equalIndifference": - equalIndifference_id.append(info[0][-3:]) - elif info[0] in subject_list and info[1] == "equalRange": - equalRange_id.append(info[0][-3:]) - - copes_equalIndifference = [] - copes_equalRange = [] - varcopes_equalIndifference = [] - varcopes_equalRange = [] - - for file in copes: - sub_id = file.split('/') - if sub_id[-2][-3:] in equalIndifference_id: - copes_equalIndifference.append(file) - elif sub_id[-2][-3:] in equalRange_id: - copes_equalRange.append(file) - - for file in varcopes: - sub_id = file.split('/') - if sub_id[-2][-3:] in equalIndifference_id: - varcopes_equalIndifference.append(file) - elif sub_id[-2][-3:] in equalRange_id: - varcopes_equalRange.append(file) - - print(len(equalRange_id)) - print(len(equalIndifference_id)) - print(len(copes_equalIndifference)) - print(len(copes_equalRange)) - - copes_global = copes_equalIndifference + copes_equalRange - varcopes_global = varcopes_equalIndifference + varcopes_equalRange - - return copes_equalIndifference, copes_equalRange, copes_global, varcopes_equalIndifference, varcopes_equalRange, varcopes_global, equalIndifference_id, equalRange_id - -def get_regs(equalRange_id, equalIndifference_id, method, subject_list): - """ - Create dictionary of regressors for group analysis. - - Parameters: - - equalRange_id: list of str, ids of subjects in equal range group - - equalIndifference_id: list of str, ids of subjects in equal indifference group - - method: one of "equalRange", "equalIndifference" or "groupComp" - - subject_list: list of str, ids of subjects for which we do the analysis - - Returns: - - regressors: dict, dictionary of regressors used to distinguish groups in FSL group analysis - """ - if method == "equalRange": - regressors = dict(group_mean = [1 for i in range(len(equalRange_id))]) - - elif method == "equalIndifference": - regressors = dict(group_mean = [1 for i in range(len(equalIndifference_id))]) - - elif method == "groupComp": - equalRange_reg = [1 for i in range(len(equalRange_id) + len(equalIndifference_id))] - equalIndifference_reg = [0 for i in range(len(equalRange_id) + len(equalIndifference_id))] - - for i, sub_id in enumerate(subject_list): - if sub_id in equalIndifference_id: - index = i - equalIndifference_reg[index] = 1 - equalRange_reg[index] = 0 - - regressors = dict(equalRange = equalRange_reg, - equalIndifference = equalIndifference_reg) - - return regressors - -def get_group_workflow(subject_list, n_sub, contrast_list, method, exp_dir, output_dir, - working_dir, result_dir, data_dir): - """ - Returns the group level of analysis workflow. - - Parameters: - - exp_dir: str, directory where raw data are stored - - result_dir: str, directory where results will be stored - - working_dir: str, name of the sub-directory for intermediate results - - output_dir: str, name of the sub-directory for final results - - subject_list: list of str, list of subject for which you want to do the preprocessing - - contrast_list: list of str, list of contrasts to analyze - - method: one of "equalRange", "equalIndifference" or "groupComp" - - n_sub: int, number of subject to include - - Returns: - - l2_analysis: Nipype WorkFlow - """ - # Infosource Node - To iterate on subject and runs - infosource_3rdlevel = Node(IdentityInterface(fields = ['contrast_id', 'exp_dir', 'result_dir', - 'output_dir', 'working_dir', 'subject_list', 'method'], - exp_dir = exp_dir, result_dir = result_dir, - output_dir = output_dir, working_dir = working_dir, - subject_list = subject_list, method = method), - name = 'infosource_3rdlevel') - infosource_3rdlevel.iterables = [('contrast_id', contrast_list)] - - # Templates to select files node - copes_file = opj(output_dir, 'l2_analysis', '_contrast_id_{contrast_id}_subject_id_*', - 'cope1.nii.gz') - - varcopes_file = opj(output_dir, 'l2_analysis', '_contrast_id_{contrast_id}_subject_id_*', - 'varcope1.nii.gz') - - participants_file = opj(exp_dir, 'participants.tsv') - - mask_file = opj(data_dir, 'NARPS-X19V', 'hypo1_unthresh.nii.gz') - - template = {'cope' : copes_file, 'varcope' : varcopes_file, 'participants' : participants_file, - 'mask':mask_file} - - # SelectFiles node - to select necessary files - selectfiles_3rdlevel = Node(SelectFiles(template, base_directory=result_dir), name = 'selectfiles_3rdlevel') - - datasink_3rdlevel = Node(DataSink(base_directory=result_dir, container=output_dir), name='datasink_3rdlevel') - - merge_copes_3rdlevel = Node(Merge(dimension = 't'), name = 'merge_copes_3rdlevel') - merge_varcopes_3rdlevel = Node(Merge(dimension = 't'), name = 'merge_varcopes_3rdlevel') - - subgroups_contrasts = Node(Function(input_names = ['copes', 'varcopes', 'subject_ids', 'participants_file'], - output_names = ['copes_equalIndifference', 'copes_equalRange', 'copes_global', - 'varcopes_equalIndifference', 'varcopes_equalRange', 'varcopes_global', - 'equalIndifference_id', 'equalRange_id'], - function = get_subgroups_contrasts), - name = 'subgroups_contrasts') - - specifymodel_3rdlevel = Node(MultipleRegressDesign(), name = 'specifymodel_3rdlevel') - - flame_3rdlevel = Node(FLAMEO(run_mode = 'flame1'), - name='flame_3rdlevel') - - 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 - - smoothest = MapNode(SmoothEstimate(), name='smoothest', iterfield = ['zstat_file']) - - cluster = MapNode(Cluster(threshold = 2.3, out_threshold_file = True, - out_pval_file = True, pthreshold=0.05), name = 'cluster', - iterfield = ['in_file', 'dlh', 'volume', 'cope_file'], synchronize = True) - - l3_analysis = Workflow(base_dir = opj(result_dir, working_dir), name = f"l3_analysis_{method}_nsub_{n_sub}") - - l3_analysis.connect([(infosource_3rdlevel, selectfiles_3rdlevel, [('contrast_id', 'contrast_id')]), - (infosource_3rdlevel, subgroups_contrasts, [('subject_list', 'subject_ids')]), - (selectfiles_3rdlevel, subgroups_contrasts, [('cope', 'copes'), ('varcope', 'varcopes'), - ('participants', 'participants_file')]), - (selectfiles_3rdlevel, flame_3rdlevel, [('mask', 'mask_file')]), - (selectfiles_3rdlevel, smoothest, [('mask', 'mask_file')]), - (subgroups_contrasts, regs, [('equalRange_id', 'equalRange_id'), - ('equalIndifference_id', 'equalIndifference_id')]), - (regs, specifymodel_3rdlevel, [('regressors', 'regressors')])]) - - - if method == 'equalIndifference' or method == 'equalRange': - specifymodel_3rdlevel.inputs.contrasts = [["group_mean", "T", ["group_mean"], [1]], ["group_mean_neg", "T", ["group_mean"], [-1]]] - - if method == 'equalIndifference': - l3_analysis.connect([(subgroups_contrasts, merge_copes_3rdlevel, - [('copes_equalIndifference', 'in_files')]), - (subgroups_contrasts, merge_varcopes_3rdlevel, - [('varcopes_equalIndifference', 'in_files')])]) - elif method == 'equalRange': - l3_analysis.connect([(subgroups_contrasts, merge_copes_3rdlevel, [('copes_equalRange', 'in_files')]), - (subgroups_contrasts, merge_varcopes_3rdlevel, [('varcopes_equalRange', 'in_files')])]) - - elif method == "groupComp": - specifymodel_3rdlevel.inputs.contrasts = [["equalRange_sup", "T", ["equalRange", "equalIndifference"], - [1, -1]]] - - l3_analysis.connect([(subgroups_contrasts, merge_copes_3rdlevel, [('copes_global', 'in_files')]), - (subgroups_contrasts, merge_varcopes_3rdlevel, [('varcopes_global', 'in_files')])]) - - l3_analysis.connect([(merge_copes_3rdlevel, flame_3rdlevel, [('merged_file', 'cope_file')]), - (merge_varcopes_3rdlevel, flame_3rdlevel, [('merged_file', 'var_cope_file')]), - (specifymodel_3rdlevel, flame_3rdlevel, [('design_mat', 'design_file'), - ('design_con', 't_con_file'), - ('design_grp', 'cov_split_file')]), - (flame_3rdlevel, cluster, [('zstats', 'in_file'), ('copes', 'cope_file')]), - (flame_3rdlevel, smoothest, [('zstats', 'zstat_file')]), - (smoothest, cluster, [('dlh', 'dlh'), ('volume', 'volume')]), - (flame_3rdlevel, datasink_3rdlevel, [('zstats', - f"l3_analysis_{method}_nsub_{n_sub}.@zstats"), - ('tstats', - f"l3_analysis_{method}_nsub_{n_sub}.@tstats")]), - (cluster, datasink_3rdlevel, [('threshold_file', f"l3_analysis_{method}_nsub_{n_sub}.@thresh"), - ('pval_file', f"l3_analysis_{method}_nsub_{n_sub}.@pval")])]) - - return l3_analysis - -def reorganize_results(result_dir, output_dir, n_sub, team_ID): - """ - Reorganize the results to analyze them. - - Parameters: - - result_dir: str, directory where results will be stored - - output_dir: str, name of the sub-directory for final results - - n_sub: float, number of subject used for the analysis - - team_ID: str, ID of the team to reorganize results - - """ - from os.path import join as opj - import os - import shutil - import gzip - import nibabel as nib - import numpy as np - - h1 = opj(result_dir, output_dir, f"l3_analysis_equalIndifference_nsub_{n_sub}", '_contrast_id_1') - h2 = opj(result_dir, output_dir, f"l3_analysis_equalRange_nsub_{n_sub}", '_contrast_id_1') - h3 = opj(result_dir, output_dir, f"l3_analysis_equalIndifference_nsub_{n_sub}", '_contrast_id_1') - h4 = opj(result_dir, output_dir, f"l3_analysis_equalRange_nsub_{n_sub}", '_contrast_id_1') - h5 = opj(result_dir, output_dir, f"l3_analysis_equalIndifference_nsub_{n_sub}", '_contrast_id_2') - h6 = opj(result_dir, output_dir, f"l3_analysis_equalRange_nsub_{n_sub}", '_contrast_id_2') - h7 = opj(result_dir, output_dir, f"l3_analysis_equalIndifference_nsub_{n_sub}", '_contrast_id_2') - h8 = opj(result_dir, output_dir, f"l3_analysis_equalRange_nsub_{n_sub}", '_contrast_id_2') - h9 = opj(result_dir, output_dir, f"l3_analysis_groupComp_nsub_{n_sub}", '_contrast_id_2') - - h = [h1, h2, h3, h4, h5, h6, h7, h8, h9] - - repro_unthresh = [opj(filename, "zstat1.nii.gz") if i in [4, 5] else opj(filename, "zstat1.nii.gz") for i, filename in enumerate(h)] - - repro_thresh = [opj(filename, '_cluster0', "zstat1_threshold.nii.gz") if i in [4, 5] else opj(filename, '_cluster0', "zstat1_threshold.nii.gz") 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.gz") - shutil.copyfile(f_in, f_out) - - for i, filename in enumerate(repro_thresh): - f_in = filename - img = nib.load(filename) - original_affine = img.affine.copy() - spm = nib.load(repro_unthresh[i]) - new_img = img.get_fdata() > 0.95 - new_img = new_img.astype(float) * spm.get_fdata() - new_spm = nib.Nifti1Image(new_img, original_affine) - nib.save(new_spm, opj(result_dir, "NARPS-reproduction", - f"team_{team_ID}_nsub_{n_sub}_hypo{i+1}_thresholded.nii.gz")) - #f_out = opj(result_dir, "final_results", f"team_{team_ID}_nsub_{n_sub}_hypo{i+1}_thresholded.nii.gz") - #shutil.copyfile(f_in, f_out) - - print(f"Results files of team {team_ID} reorganized.") + onsets['gain'].append(float(info[0])) + durations['gain'].append(4.0) + amplitudes['gain'].append(float(info[2])) + onsets['loss'].append(float(info[0])) + durations['loss'].append(4.0) + amplitudes['loss'].append(float(info[3])) + onsets['trial'].append(float(info[0])) + durations['trial'].append(4.0) + amplitudes['trial'].append(1.0) + + amplitudes['gain'] = amplitudes['gain'] - mean(amplitudes['gain']) + amplitudes['loss'] = amplitudes['loss'] - mean(amplitudes['loss']) + + return [ + Bunch( + conditions = condition_names, + onsets = [onsets[k] for k in condition_names], + durations = [durations[k] for k in condition_names], + amplitudes = [amplitudes[k] for k in condition_names], + regressor_names = None, + regressors = None) + ] + + def get_confounds_file(filepath, subject_id, run_id, working_dir): + """ + 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 + - working_dir: str, name of the directory for intermediate results + + Return : + - confounds_file : paths to new files containing only desired confounds. + """ + from os import makedirs + from os.path import join + + from pandas import read_csv, DataFrame + from numpy import array, transpose + + data_frame = read_csv(filepath, sep = '\t', header=0) + temp_list = array([ + data_frame['X'], data_frame['Y'], data_frame['Z'], + data_frame['RotX'], data_frame['RotY'], data_frame['RotZ']]) + retained_confounds = DataFrame(transpose(temp_list)) + + 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') 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_analysis : nipype.WorkFlow + """ + # IdentityInterface Node - To 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 - to select 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'), + 'event' : join('sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-{run_id}_events.tsv'), + 'param' : join('derivatives', 'fmriprep', 'sub-{subject_id}', 'func', + 'sub-{subject_id}_task-MGT_run-{run_id}_bold_confounds.tsv') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory = self.directories.dataset_dir + + # 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 + + # BET Node - Skullstripping data + skull_stripping_func = Node(BET(), name = 'skull_stripping_func') + skull_stripping_func.inputs.frac = 0.1 + skull_stripping_func.inputs.functional = True + skull_stripping_func.inputs.mask = True + + # IsotropicSmooth Node - Smoothing data + smoothing_func = Node(IsotropicSmooth(), name = 'smoothing_func') + smoothing_func.inputs.fwhm = self.fwhm + + # 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') + + # SpecifyModel Node - Generate run level model + specify_model = Node(SpecifyModel(), name = 'specify_model') + specify_model.inputs.high_pass_filter_cutoff = 60 + specify_model.inputs.input_units = 'secs' + specify_model.inputs.time_repetition = TaskInformation()['RepetitionTime'] + + # Function Node get_confounds_file - Get files with movement confounds + confounds = Node(Function( + function = self.get_confounds_file, + input_names = ['filepath', 'subject_id', 'run_id', 'working_dir'], + output_names = ['confounds_file']), + name = 'confounds') + confounds.inputs.working_dir = self.directories.working_dir + + # Level1Design Node - Generate files for run level computation + model_design = Node(Level1Design(), name = 'model_design') + model_design.inputs.bases = {'dgamma':{'derivs' : False}} + model_design.inputs.interscan_interval = TaskInformation()['RepetitionTime'] + model_design.inputs.model_serial_correlations = True + model_design.inputs.contrasts = self.run_level_contrasts + + # FEATModel Node - Generate run level model + model_generation = Node(FEATModel(), name = 'model_generation') + + # FILMGLS Node - Estimate first level model + model_estimate = Node(FILMGLS(), name='model_estimate') + + # Create run level analysis workflow and connect its nodes + run_level_analysis = Workflow( + base_dir = self.directories.working_dir, + name = 'run_level_analysis' + ) + run_level_analysis.connect([ + (information_source, select_files, [ + ('subject_id', 'subject_id'), + ('run_id', 'run_id')]), + (select_files, subject_information, [('event', 'event_file')]), + (select_files, confounds, [('param', 'filepath')]), + (information_source, confounds, [ + ('subject_id', 'subject_id'), + ('run_id', 'run_id')]), + (select_files, skull_stripping_func, [('func', 'in_file')]), + (skull_stripping_func, smoothing_func, [('out_file', 'in_file')]), + (confounds, specify_model, [('confounds_file', 'realignment_parameters')]), + (smoothing_func, specify_model, [('out_file', 'functional_runs')]), + (subject_information, specify_model, [('subject_info', 'subject_info')]), + (specify_model, model_design, [('session_info', 'session_info')]), + (model_design, model_generation, [ + ('ev_files', 'ev_files'), + ('fsf_files', 'fsf_file')]), + (smoothing_func, model_estimate, [('out_file', 'in_file')]), + (model_generation, model_estimate, [ + ('con_file', 'tcon_file'), + ('design_file', 'design_file')]), + (model_estimate, data_sink, [('results_dir', 'run_level_analysis.@results')]), + (model_generation, data_sink, [ + ('design_file', 'run_level_analysis.@design_file'), + ('design_image', 'run_level_analysis.@design_img')]), + (skull_stripping_func, data_sink, [('mask_file', 'run_level_analysis.@skullstriped')]) + ]) + + # Remove large files, if requested + if Configuration()['pipelines']['remove_unused_data']: + + # Remove Node - Remove skullstriped func files once they are no longer needed + remove_skullstrip = Node( + InterfaceFactory.create('remove_parent_directory'), + name = 'remove_skullstrip') + + # Remove Node - Remove smoothed files once they are no longer needed + remove_smooth = Node( + InterfaceFactory.create('remove_parent_directory'), + name = 'remove_smooth') + + # Add connections + run_level_analysis.connect([ + (data_sink, remove_skullstrip, [('out_file', '_')]), + (skull_stripping_func, remove_skullstrip, [('out_file', 'file_name')]), + (data_sink, remove_smooth, [('out_file', '_')]), + (smoothing_func, remove_smooth, [('out_file', 'file_name')]) + ]) + + return run_level_analysis + + 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, + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'run_level_analysis', '_run_id_{run_id}_subject_id_{subject_id}', + 'sub-{subject_id}_task-MGT_run-{run_id}_bold_space-MNI152NLin2009cAsym_preproc_brain_mask.nii.gz' + ) + return_list = [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + 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_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets for template in templates] + + return return_list + + def get_subject_level_analysis(self): + """ + Create the subject level analysis workflow. + + Returns: + - subject_level_analysis : nipype.WorkFlow + """ + + # 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(self.directories.output_dir, + 'run_level_analysis', '_run_id_*_subject_id_{subject_id}', + 'sub-{subject_id}_task-MGT_run-*_bold_space-MNI152NLin2009cAsym_preproc_brain_mask.nii.gz') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory= self.directories.results_dir + + # 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 + + # L2Model Node - Generate subject specific second level model + generate_model = Node(L2Model(), name = 'generate_model') + generate_model.inputs.num_copes = len(self.run_list) + + # Merge Node - Merge copes files for each subject + merge_copes = Node(Merge(), name = 'merge_copes') + merge_copes.inputs.dimension = 't' + + # Merge Node - Merge varcopes files for each subject + merge_varcopes = Node(Merge(), name = 'merge_varcopes') + merge_varcopes.inputs.dimension = 't' + + # 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 + + # 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) + + # FLAMEO Node - Estimate model + estimate_model = Node(FLAMEO(), name = 'estimate_model') + estimate_model.inputs.run_mode = 'flame1' + + # Second level (single-subject, mean of all four scans) analysis workflow. + subject_level_analysis = Workflow( + base_dir = self.directories.working_dir, + name = 'subject_level_analysis') + subject_level_analysis.connect([ + (information_source, select_files, [ + ('subject_id', 'subject_id'), + ('contrast_id', 'contrast_id')]), + (select_files, merge_copes, [('cope', 'in_files')]), + (select_files, merge_varcopes, [('varcope', 'in_files')]), + (select_files, split_masks, [('masks', 'inlist')]), + (split_masks, mask_intersection, [('out1', 'in_file')]), + (split_masks, mask_intersection, [('out2', 'operand_files')]), + (mask_intersection, estimate_model, [('out_file', 'mask_file')]), + (merge_copes, estimate_model, [('merged_file', 'cope_file')]), + (merge_varcopes, estimate_model, [('merged_file', 'var_cope_file')]), + (generate_model, estimate_model, [ + ('design_mat', 'design_file'), + ('design_con', 't_con_file'), + ('design_grp', 'cov_split_file')]), + (mask_intersection, data_sink, [('out_file', 'subject_level_analysis.@mask')]), + (estimate_model, data_sink, [ + ('zstats', 'subject_level_analysis.@stats'), + ('tstats', 'subject_level_analysis.@tstats'), + ('copes', 'subject_level_analysis.@copes'), + ('var_copes', 'subject_level_analysis.@varcopes')])]) + + return subject_level_analysis + + def get_subject_level_outputs(self): + """ Return the names of the files the subject level analysis is supposed to generate. """ + + parameters = { + 'contrast_id' : self.contrast_list, + 'subject_id' : self.subject_list, + 'file' : ['cope1.nii.gz', 'tstat1.nii.gz', 'varcope1.nii.gz', 'zstat1.nii.gz'] + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'subject_level_analysis', '_contrast_id_{contrast_id}_subject_id_{subject_id}','{file}' + ) + return_list = [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + parameters = { + 'contrast_id' : self.contrast_list, + 'subject_id' : self.subject_list, + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'subject_level_analysis', '_contrast_id_{contrast_id}_subject_id_{subject_id}', + 'sub-{subject_id}_task-MGT_run-01_bold_space-MNI152NLin2009cAsym_preproc_brain_mask_maths.nii.gz' + ) + return_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + return return_list + + def get_one_sample_t_test_regressors(subject_list: list) -> dict: + """ + Create dictionary of regressors for one sample t-test group analysis. + + Parameters: + - subject_list: ids of subject in the group for which to do the analysis + + Returns: + - dict containing named lists of regressors. + """ + + return dict(group_mean = [1 for _ in subject_list]) + + def get_two_sample_t_test_regressors( + equal_range_ids: list, + equal_indifference_ids: list, + subject_list: list, + ) -> dict: + """ + Create dictionary of regressors for two sample t-test group analysis. + + Parameters: + - equal_range_ids: ids of subjects in equal range group + - equal_indifference_ids: ids of subjects in equal indifference group + - subject_list: ids of subject for which to do the analysis + + Returns: + - regressors, dict: containing named lists of regressors. + - groups, list: group identifiers to distinguish groups in FSL analysis. + """ + + # Create 2 lists containing n_sub values which are + # * 1 if the participant is on the group + # * 0 otherwise + equal_range_regressors = [1 if i in equal_range_ids else 0 for i in subject_list] + equal_indifference_regressors = [ + 1 if i in equal_indifference_ids else 0 for i in subject_list + ] + + # Create regressors output : a dict with the two list + regressors = dict( + equalRange = equal_range_regressors, + equalIndifference = equal_indifference_regressors + ) + + # Create groups outputs : a list with 1 for equalRange subjects and 2 for equalIndifference + groups = [1 if i == 1 else 2 for i in equal_range_regressors] + + return regressors, groups + + def get_group_level_analysis(self): + """ + Return all workflows for the group level analysis. + + Returns; + - a list of nipype.WorkFlow + """ + + methods = ['equalRange', 'equalIndifference', 'groupComp'] + return [self.get_group_level_analysis_sub_workflow(method) for method in methods] + + def get_group_level_analysis_sub_workflow(self, method): + """ + Return a workflow for the group level analysis. + + Parameters: + - method: one of 'equalRange', 'equalIndifference' or 'groupComp' + + Returns: + - group_level_analysis: nipype.WorkFlow + """ + # Infosource Node - iterate over the contrasts generated by the subject level analysis + information_source = Node(IdentityInterface( + fields = ['contrast_id']), + name = 'information_source') + information_source.iterables = [('contrast_id', self.contrast_list)] + + # SelectFiles Node - select necessary files + templates = { + 'cope' : join(self.directories.output_dir, + 'subject_level_analysis', '_contrast_id_{contrast_id}_subject_id_*', + 'cope1.nii.gz'), + 'varcope' : join(self.directories.output_dir, + 'subject_level_analysis', '_contrast_id_{contrast_id}_subject_id_*', + 'varcope1.nii.gz'), + 'masks': join(self.directories.output_dir, + 'subject_level_analysis', '_contrast_id_1_subject_id_*', + 'sub-*_task-MGT_run-*_bold_space-MNI152NLin2009cAsym_preproc_brain_mask_maths.nii.gz') + } + select_files = Node(SelectFiles(templates), name = 'select_files') + select_files.inputs.base_directory = self.directories.results_dir + + # Datasink Node - save important files + data_sink = Node(DataSink(), name = 'data_sink') + data_sink.inputs.base_directory = self.directories.output_dir + + # 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' + ) + + # 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' + ) + + # Merge Node - Merge cope files + merge_copes = Node(Merge(), name = 'merge_copes') + merge_copes.inputs.dimension = 't' + + # Merge Node - Merge cope files + merge_varcopes = Node(Merge(), name = 'merge_varcopes') + merge_varcopes.inputs.dimension = 't' + + # 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 + + # 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) + + # 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' + + # SmoothEstimate Node - Estimates the smoothness of zstat files + smoothness_estimate = MapNode(SmoothEstimate(), + name='smoothness_estimate', + iterfield = ['zstat_file']) + + # Cluster Node - Perform clustering on statistical output + cluster = MapNode(Cluster(), + name = 'cluster', + iterfield = ['in_file', 'dlh', 'volume', 'cope_file'], + synchronize = True) + cluster.inputs.threshold = 2.3 + cluster.inputs.out_threshold_file = True + cluster.inputs.out_pval_file = True + cluster.inputs.pthreshold = 0.05 + + # Compute the number of participants used to do the analysis + nb_subjects = len(self.subject_list) + + # Declare the workflow + 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_copes, [('cope', 'input_str')]), + (select_files, get_varcopes, [('varcope', 'input_str')]), + (get_copes, merge_copes, [(('out_list', clean_list), 'in_files')]), + (get_varcopes, merge_varcopes,[(('out_list', clean_list), 'in_files')]), + (select_files, split_masks, [('masks', 'inlist')]), + (split_masks, mask_intersection, [('out1', 'in_file')]), + (split_masks, mask_intersection, [('out2', 'operand_files')]), + (mask_intersection, estimate_model, [('out_file', 'mask_file')]), + (mask_intersection, smoothness_estimate, [('out_file', 'mask_file')]), + (merge_copes, estimate_model, [('merged_file', 'cope_file')]), + (merge_varcopes, estimate_model, [('merged_file', 'var_cope_file')]), + (specify_model, estimate_model, [ + ('design_mat', 'design_file'), + ('design_con', 't_con_file'), + ('design_grp', 'cov_split_file') + ]), + (estimate_model, cluster, [('zstats', 'in_file'), ('copes', 'cope_file')]), + (estimate_model, smoothness_estimate, [('zstats', 'zstat_file')]), + (smoothness_estimate, cluster, [('dlh', 'dlh'), ('volume', 'volume')]), + (estimate_model, data_sink, [ + ('zstats', f'group_level_analysis_{method}_nsub_{nb_subjects}.@zstats'), + ('tstats', f'group_level_analysis_{method}_nsub_{nb_subjects}.@tstats') + ]), + (cluster, data_sink, [ + ('threshold_file', f'group_level_analysis_{method}_nsub_{nb_subjects}.@thresh'), + ('pval_file', f'group_level_analysis_{method}_nsub_{nb_subjects}.@pval')]) + ]) + + if method in ('equalIndifference', 'equalRange'): + # Setup a one sample t-test + specify_model.inputs.contrasts = [ + ['group_mean', 'T', ['group_mean'], [1]], + ['group_mean_neg', 'T', ['group_mean'], [-1]] + ] + + # Function Node get_group_subjects - Get subjects in the group and in the subject_list + get_group_subjects = Node(Function( + function = list_intersection, + input_names = ['list_1', 'list_2'], + output_names = ['out_list'] + ), + name = 'get_group_subjects' + ) + get_group_subjects.inputs.list_1 = get_group(method) + get_group_subjects.inputs.list_2 = self.subject_list + + # Function Node get_one_sample_t_test_regressors + # Get regressors in the equalRange and equalIndifference method case + regressors_one_sample = Node( + Function( + function = self.get_one_sample_t_test_regressors, + input_names = ['subject_list'], + output_names = ['regressors'] + ), + name = 'regressors_one_sample', + ) + + # Add missing connections + group_level_analysis.connect([ + (get_group_subjects, get_copes, [('out_list', 'elements')]), + (get_group_subjects, get_varcopes, [('out_list', 'elements')]), + (get_group_subjects, regressors_one_sample, [('out_list', 'subject_list')]), + (regressors_one_sample, specify_model, [('regressors', 'regressors')]) + ]) + + elif method == 'groupComp': + + # Select copes and varcopes corresponding to the selected subjects + # Indeed the SelectFiles node asks for all (*) subjects available + get_copes.inputs.elements = self.subject_list + get_varcopes.inputs.elements = self.subject_list + + # Setup a two sample t-test + specify_model.inputs.contrasts = [ + ['equalRange_sup', 'T', ['equalRange', 'equalIndifference'], [1, -1]] + ] + + # Function Node get_equal_range_subjects + # Get subjects in the equalRange group and in the subject_list + get_equal_range_subjects = Node(Function( + function = list_intersection, + input_names = ['list_1', 'list_2'], + output_names = ['out_list'] + ), + name = 'get_equal_range_subjects' + ) + get_equal_range_subjects.inputs.list_1 = get_group('equalRange') + get_equal_range_subjects.inputs.list_2 = self.subject_list + + # Function Node get_equal_indifference_subjects + # Get subjects in the equalIndifference group and in the subject_list + get_equal_indifference_subjects = Node(Function( + function = list_intersection, + input_names = ['list_1', 'list_2'], + output_names = ['out_list'] + ), + name = 'get_equal_indifference_subjects' + ) + get_equal_indifference_subjects.inputs.list_1 = get_group('equalIndifference') + get_equal_indifference_subjects.inputs.list_2 = self.subject_list + + # Function Node get_two_sample_t_test_regressors + # Get regressors in the groupComp method case + regressors_two_sample = Node( + Function( + function = self.get_two_sample_t_test_regressors, + input_names = [ + 'equal_range_ids', + 'equal_indifference_ids', + 'subject_list', + ], + output_names = ['regressors', 'groups'] + ), + name = 'regressors_two_sample', + ) + regressors_two_sample.inputs.subject_list = self.subject_list + + # Add missing connections + group_level_analysis.connect([ + (get_equal_range_subjects, regressors_two_sample, [ + ('out_list', 'equal_range_ids') + ]), + (get_equal_indifference_subjects, regressors_two_sample, [ + ('out_list', 'equal_indifference_ids') + ]), + (regressors_two_sample, specify_model, [ + ('regressors', 'regressors'), + ('groups', 'groups')]) + ]) + + return group_level_analysis + + def get_group_level_outputs(self): + """ Return all names for the files the group level analysis is supposed to generate. """ + + # Handle equalRange and equalIndifference + parameters = { + 'contrast_id': self.contrast_list, + 'method': ['equalRange', 'equalIndifference'], + 'file': [ + '_cluster0/zstat1_pval.nii.gz ', + '_cluster0/zstat1_threshold.nii.gz', + '_cluster1/zstat2_pval.nii.gz', + '_cluster1/zstat2_threshold.nii.gz', + 'tstat1.nii.gz', + 'tstat2.nii.gz', + 'zstat1.nii.gz', + 'zstat2.nii.gz' + ] + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + 'group_level_analysis_{method}_nsub_'+f'{len(self.subject_list)}', + '_contrast_id_{contrast_id}', + '{file}' + ) + return_list = [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + # Handle groupComp + parameters = { + 'contrast_id': self.contrast_list, + 'file': [ + '_cluster0/zstat1_pval.nii.gz ', + '_cluster0/zstat1_threshold.nii.gz', + 'tstat1.nii.gz', + 'zstat1.nii.gz' + ] + } + parameter_sets = product(*parameters.values()) + template = join( + self.directories.output_dir, + f'group_level_analysis_groupComp_nsub_{len(self.subject_list)}', + '_contrast_id_{contrast_id}', + '{file}' + ) + return_list += [template.format(**dict(zip(parameters.keys(), parameter_values)))\ + for parameter_values in parameter_sets] + + return return_list + + def get_hypotheses_outputs(self): + """ Return all hypotheses output file names. """ + + nb_sub = len(self.subject_list) + files = [ + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_1', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalIndifference_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_equalRange_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz'), + join(f'group_level_analysis_groupComp_nsub_{nb_sub}', + '_contrast_id_2', '_cluster0', 'zstat1_threshold.nii.gz'), + join(f'group_level_analysis_groupComp_nsub_{nb_sub}', + '_contrast_id_2', 'zstat1.nii.gz') + ] + return [join(self.directories.output_dir, f) for f in files] diff --git a/tests/conftest.py b/tests/conftest.py index 680dea31..5f78f36f 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -10,6 +10,7 @@ from os.path import join, isfile from shutil import rmtree +from numpy import isclose from pytest import helpers from pathvalidate import is_valid_filepath from numpy import isclose diff --git a/tests/core/test_common.py b/tests/core/test_common.py index 3bfb69f8..1c362023 100644 --- a/tests/core/test_common.py +++ b/tests/core/test_common.py @@ -59,7 +59,7 @@ def test_remove_file(remove_test_dir): # Check file is removed assert not exists(test_file_path) - + @staticmethod @mark.unit_test def test_remove_directory(remove_test_dir): @@ -396,7 +396,7 @@ def test_node_list_to_file_1(): # Check file was created assert exists(out_file) - + # Check file was created with open(out_file, 'r', encoding = 'utf-8') as file: for list_element, file_element in zip(out_list, file.read().split('\n')): @@ -429,9 +429,8 @@ def test_node_list_to_file_2(): # Check file was created assert exists(out_file) - + # Check file was created with open(out_file, 'r', encoding = 'utf-8') as file: for list_element, file_element in zip(out_list, file.read().split('\n')): assert list_element == file_element - diff --git a/tests/data/test_results.py b/tests/data/test_results.py index 2465eb7d..d1a7a93c 100644 --- a/tests/data/test_results.py +++ b/tests/data/test_results.py @@ -12,8 +12,7 @@ """ from os.path import isdir, join -from shutil import rmtree, move, copytree -from time import sleep +from shutil import rmtree, copytree from checksumdir import dirhash from pytest import mark @@ -74,7 +73,8 @@ def fake_get_uid(_): # Mock the results path results_directory = Configuration()['directories']['narps_results'] - Configuration()['directories']['narps_results'] = Configuration()['directories']['test_runs'] + Configuration()[ + 'directories']['narps_results'] = Configuration()['directories']['test_runs'] # Init & download the collection collection = ResultsCollection('2T6S') diff --git a/tests/data/test_task.py b/tests/data/test_task.py index 8b6860dd..229a1ddc 100644 --- a/tests/data/test_task.py +++ b/tests/data/test_task.py @@ -15,7 +15,7 @@ from pytest import mark, fixture from narps_open.utils.configuration import Configuration -import narps_open.data.task as task +from narps_open.data import task @fixture(scope='function', autouse=True) def mock_task_data(mocker): @@ -50,7 +50,7 @@ def test_singleton(): @mark.unit_test def test_derived(): """ Test the derived values of a TaskInformation object """ - + task_info = task.TaskInformation() assert task_info['NumberOfSlices'] == 6 assert task_info['AcquisitionTime'] == 1 / 6 diff --git a/tests/pipelines/test_team_08MQ.py b/tests/pipelines/test_team_08MQ.py index b962557f..09fdc66f 100644 --- a/tests/pipelines/test_team_08MQ.py +++ b/tests/pipelines/test_team_08MQ.py @@ -57,19 +57,11 @@ def test_outputs(): pipeline = PipelineTeam08MQ() # 1 - 1 subject outputs pipeline.subject_list = ['001'] - assert len(pipeline.get_preprocessing_outputs()) == 4*4 - assert len(pipeline.get_run_level_outputs()) == 8+4*3*4 - assert len(pipeline.get_subject_level_outputs()) == 4*3 - assert len(pipeline.get_group_level_outputs()) == 0 - assert len(pipeline.get_hypotheses_outputs()) == 18 + helpers.test_pipeline_outputs(pipeline, [4*4, 8+4*3*4, 4*3, 8*2*3 + 4*3, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - assert len(pipeline.get_preprocessing_outputs()) == 4*4*4 - assert len(pipeline.get_run_level_outputs()) == (8+4*3*4)*4 - assert len(pipeline.get_subject_level_outputs()) == 4*3*4 - assert len(pipeline.get_group_level_outputs()) == 0 - assert len(pipeline.get_hypotheses_outputs()) == 18 + helpers.test_pipeline_outputs(pipeline, [4*4*4, (8+4*3*4)*4, 4*3*4, 8*2*3 + 4*3, 18]) @staticmethod @mark.unit_test @@ -130,10 +122,10 @@ def test_two_sample_t_test_regressors(): ['002', '004'], # equalIndifference group ['001', '002', '003', '004'] # all subjects ) - assert regressors == dict( - equalRange = [1, 0, 1, 0], - equalIndifference = [0, 1, 0, 1] - ) + assert regressors == { + 'equalRange' : [1, 0, 1, 0], + 'equalIndifference' : [0, 1, 0, 1] + } assert groups == [1, 2, 1, 2] @staticmethod diff --git a/tests/pipelines/test_team_2T6S.py b/tests/pipelines/test_team_2T6S.py index 43c41773..fbe5a9ad 100644 --- a/tests/pipelines/test_team_2T6S.py +++ b/tests/pipelines/test_team_2T6S.py @@ -47,19 +47,11 @@ def test_outputs(): pipeline = PipelineTeam2T6S() # 1 - 1 subject outputs pipeline.subject_list = ['001'] - assert len(pipeline.get_preprocessing_outputs()) == 0 - assert len(pipeline.get_run_level_outputs()) == 0 - assert len(pipeline.get_subject_level_outputs()) == 7 - assert len(pipeline.get_group_level_outputs()) == 63 - assert len(pipeline.get_hypotheses_outputs()) == 18 + helpers.test_pipeline_outputs(pipeline, [0, 0, 7, 63, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - assert len(pipeline.get_preprocessing_outputs()) == 0 - assert len(pipeline.get_run_level_outputs()) == 0 - assert len(pipeline.get_subject_level_outputs()) == 28 - assert len(pipeline.get_group_level_outputs()) == 63 - assert len(pipeline.get_hypotheses_outputs()) == 18 + helpers.test_pipeline_outputs(pipeline, [0, 0, 28, 63, 18]) @staticmethod @mark.pipeline_test diff --git a/tests/pipelines/test_team_C88N.py b/tests/pipelines/test_team_C88N.py index e6bb70fd..7fbdbe3e 100644 --- a/tests/pipelines/test_team_C88N.py +++ b/tests/pipelines/test_team_C88N.py @@ -52,19 +52,11 @@ def test_outputs(): pipeline = PipelineTeamC88N() # 1 - 1 subject outputs pipeline.subject_list = ['001'] - assert len(pipeline.get_preprocessing_outputs()) == 0 - assert len(pipeline.get_run_level_outputs()) == 0 - assert len(pipeline.get_subject_level_outputs()) == 8 - assert len(pipeline.get_group_level_outputs()) == 53 - assert len(pipeline.get_hypotheses_outputs()) == 18 + helpers.test_pipeline_outputs(pipeline, [0, 0, 8, 53, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - assert len(pipeline.get_preprocessing_outputs()) == 0 - assert len(pipeline.get_run_level_outputs()) == 0 - assert len(pipeline.get_subject_level_outputs()) == 32 - assert len(pipeline.get_group_level_outputs()) == 53 - assert len(pipeline.get_hypotheses_outputs()) == 18 + helpers.test_pipeline_outputs(pipeline, [0, 0, 32, 53, 18]) @staticmethod @mark.unit_test @@ -72,7 +64,8 @@ def test_subject_information(): """ Test the get_subject_information method """ # Test with 'gain' - test_event_file = join(Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') + test_event_file = join( + Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') information = PipelineTeamC88N.get_subject_information( [test_event_file, test_event_file], 'gain' @@ -103,7 +96,8 @@ def test_subject_information(): assert isclose(reference_array, test_array).all() # Test with 'loss' - test_event_file = join(Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') + test_event_file = join( + Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') information = PipelineTeamC88N.get_subject_information( [test_event_file, test_event_file], 'loss' diff --git a/tests/pipelines/test_team_J7F9.py b/tests/pipelines/test_team_J7F9.py index 706c0b56..da37e990 100644 --- a/tests/pipelines/test_team_J7F9.py +++ b/tests/pipelines/test_team_J7F9.py @@ -16,7 +16,6 @@ from filecmp import cmp from pytest import helpers, mark, fixture -from numpy import isclose from nipype import Workflow from nipype.interfaces.base import Bunch @@ -34,14 +33,6 @@ def remove_test_dir(): yield # test runs here rmtree(TEMPORARY_DIR, ignore_errors = True) -def compare_float_2d_arrays(array_1, array_2): - """ Assert array_1 and array_2 are close enough """ - - assert len(array_1) == len(array_2) - for reference_array, test_array in zip(array_1, array_2): - assert len(reference_array) == len(test_array) - assert isclose(reference_array, test_array).all() - class TestPipelinesTeamJ7F9: """ A class that contains all the unit tests for the PipelineTeamJ7F9 class.""" @@ -73,19 +64,11 @@ def test_outputs(): pipeline = PipelineTeamJ7F9() # 1 - 1 subject outputs pipeline.subject_list = ['001'] - assert len(pipeline.get_preprocessing_outputs()) == 0 - assert len(pipeline.get_run_level_outputs()) == 0 - assert len(pipeline.get_subject_level_outputs()) == 7 - assert len(pipeline.get_group_level_outputs()) == 63 - assert len(pipeline.get_hypotheses_outputs()) == 18 + helpers.test_pipeline_outputs(pipeline, [0, 0, 7, 63, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - assert len(pipeline.get_preprocessing_outputs()) == 0 - assert len(pipeline.get_run_level_outputs()) == 0 - assert len(pipeline.get_subject_level_outputs()) == 28 - assert len(pipeline.get_group_level_outputs()) == 63 - assert len(pipeline.get_hypotheses_outputs()) == 18 + helpers.test_pipeline_outputs(pipeline, [0, 0, 28, 63, 18]) @staticmethod @mark.unit_test @@ -106,80 +89,102 @@ def test_subject_information(): bunch = info_missed[0] assert isinstance(bunch, Bunch) assert bunch.conditions == ['trial', 'missed'] - compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 19.535, 27.535, 36.435], [19.535]]) - compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0]]) - assert bunch.amplitudes == None - assert bunch.tmod == None + helpers.compare_float_2d_arrays(bunch.onsets, [ + [4.071, 11.834, 19.535, 27.535, 36.435], [19.535]]) + helpers.compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.0, 0.0], [0.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] - compare_float_2d_arrays(bunch.pmod[0].param, [[-8.4, 11.6, 15.6, -12.4, -6.4], [-8.2, -0.2, 4.8, 0.8, 2.8]]) - assert bunch.regressor_names == None - assert bunch.regressors == None + helpers.compare_float_2d_arrays(bunch.pmod[0].param, [ + [-8.4, 11.6, 15.6, -12.4, -6.4], [-8.2, -0.2, 4.8, 0.8, 2.8]]) + assert bunch.regressor_names is None + assert bunch.regressors is None bunch = info_missed[1] assert isinstance(bunch, Bunch) assert bunch.conditions == ['trial', 'missed'] - compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 19.535, 27.535, 36.435], [19.535]]) - compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0]]) - assert bunch.amplitudes == None - assert bunch.tmod == None + helpers.compare_float_2d_arrays(bunch.onsets, [ + [4.071, 11.834, 19.535, 27.535, 36.435], [19.535]]) + helpers.compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.0, 0.0], [0.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] - compare_float_2d_arrays(bunch.pmod[0].param, [[-8.4, 11.6, 15.6, -12.4, -6.4], [-8.2, -0.2, 4.8, 0.8, 2.8]]) - assert bunch.regressor_names == None - assert bunch.regressors == None + helpers.compare_float_2d_arrays(bunch.pmod[0].param, [ + [-8.4, 11.6, 15.6, -12.4, -6.4], [-8.2, -0.2, 4.8, 0.8, 2.8]]) + assert bunch.regressor_names is None + assert bunch.regressors is None bunch = info_ok[0] assert isinstance(bunch, Bunch) assert bunch.conditions == ['trial'] - compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 27.535, 36.435]]) - compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.0]]) - assert bunch.amplitudes == None - assert bunch.tmod == None + helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 27.535, 36.435]]) + helpers.compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.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] - compare_float_2d_arrays(bunch.pmod[0].param, [[-4.5, 15.5, -8.5, -2.5], [-7.0, 1.0, 2.0, 4.0]]) - assert bunch.regressor_names == None - assert bunch.regressors == None + helpers.compare_float_2d_arrays(bunch.pmod[0].param, [ + [-4.5, 15.5, -8.5, -2.5], [-7.0, 1.0, 2.0, 4.0]]) + assert bunch.regressor_names is None + assert bunch.regressors is None bunch = info_ok[1] assert isinstance(bunch, Bunch) assert bunch.conditions == ['trial'] - compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 27.535, 36.435]]) - compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.0]]) - assert bunch.amplitudes == None - assert bunch.tmod == None + helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 27.535, 36.435]]) + helpers.compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.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] - compare_float_2d_arrays(bunch.pmod[0].param, [[-4.5, 15.5, -8.5, -2.5], [-7.0, 1.0, 2.0, 4.0]]) - assert bunch.regressor_names == None - assert bunch.regressors == None - + helpers.compare_float_2d_arrays(bunch.pmod[0].param, [ + [-4.5, 15.5, -8.5, -2.5], [-7.0, 1.0, 2.0, 4.0]]) + assert bunch.regressor_names is None + assert bunch.regressors is None + bunch = info_half[0] assert isinstance(bunch, Bunch) assert bunch.conditions == ['trial', 'missed'] - compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 27.535, 36.435], []]) - compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.0], []]) - assert bunch.amplitudes == None - assert bunch.tmod == None + helpers.compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 27.535, 36.435], []]) + helpers.compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.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] - compare_float_2d_arrays(bunch.pmod[0].param, [[-6.666666666666668, 13.333333333333332, -10.666666666666668, -4.666666666666668], [-7.666666666666666, 0.3333333333333339, 1.333333333333334, 3.333333333333334]]) - assert bunch.regressor_names == None - assert bunch.regressors == None + helpers.compare_float_2d_arrays( + bunch.pmod[0].param, + [[-6.666666666666668, 13.333333333333332, -10.666666666666668, -4.666666666666668], + [-7.666666666666666, 0.3333333333333339, 1.333333333333334, 3.333333333333334]]) + assert bunch.regressor_names is None + assert bunch.regressors is None bunch = info_half[1] assert isinstance(bunch, Bunch) assert bunch.conditions == ['trial', 'missed'] - compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 19.535, 27.535, 36.435], [19.535]]) - compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.0, 0.0], [0.0]]) - assert bunch.amplitudes == None - assert bunch.tmod == None + helpers.compare_float_2d_arrays(bunch.onsets, [ + [4.071, 11.834, 19.535, 27.535, 36.435], [19.535]]) + helpers.compare_float_2d_arrays(bunch.durations, [[0.0, 0.0, 0.0, 0.0, 0.0], [0.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] - compare_float_2d_arrays(bunch.pmod[0].param, [[-6.666666666666668, 13.333333333333332, 17.333333333333332, -10.666666666666668, -4.666666666666668], [-7.666666666666666, 0.3333333333333339, 5.333333333333334, 1.333333333333334, 3.333333333333334]]) - assert bunch.regressor_names == None - assert bunch.regressors == None + helpers.compare_float_2d_arrays( + bunch.pmod[0].param, + [[ + -6.666666666666668, + 13.333333333333332, + 17.333333333333332, + -10.666666666666668, + -4.666666666666668], + [-7.666666666666666, + 0.3333333333333339, + 5.333333333333334, + 1.333333333333334, + 3.333333333333334]]) + assert bunch.regressor_names is None + assert bunch.regressors is None @staticmethod @mark.unit_test diff --git a/tests/pipelines/test_team_Q6O0.py b/tests/pipelines/test_team_Q6O0.py index 639f609e..6fc7513c 100644 --- a/tests/pipelines/test_team_Q6O0.py +++ b/tests/pipelines/test_team_Q6O0.py @@ -47,19 +47,11 @@ def test_outputs(): pipeline = PipelineTeamQ6O0() # 1 - 1 subject outputs pipeline.subject_list = ['001'] - assert len(pipeline.get_preprocessing_outputs()) == 0 - assert len(pipeline.get_run_level_outputs()) == 0 - assert len(pipeline.get_subject_level_outputs()) == 6 - assert len(pipeline.get_group_level_outputs()) == 37 - assert len(pipeline.get_hypotheses_outputs()) == 18 + helpers.test_pipeline_outputs(pipeline, [0, 0, 6, 37, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - assert len(pipeline.get_preprocessing_outputs()) == 0 - assert len(pipeline.get_run_level_outputs()) == 0 - assert len(pipeline.get_subject_level_outputs()) == 24 - assert len(pipeline.get_group_level_outputs()) == 37 - assert len(pipeline.get_hypotheses_outputs()) == 18 + helpers.test_pipeline_outputs(pipeline, [0, 0, 24, 37, 18]) @staticmethod @mark.pipeline_test diff --git a/tests/pipelines/test_team_T54A.py b/tests/pipelines/test_team_T54A.py index 05ecc003..2199cf40 100644 --- a/tests/pipelines/test_team_T54A.py +++ b/tests/pipelines/test_team_T54A.py @@ -15,7 +15,6 @@ from shutil import rmtree from pytest import helpers, mark, fixture -from numpy import isclose from nipype import Workflow from nipype.interfaces.base import Bunch @@ -33,14 +32,6 @@ def remove_test_dir(): yield # test runs here rmtree(TEMPORARY_DIR, ignore_errors = True) -def compare_float_2d_arrays(array_1, array_2): - """ Assert array_1 and array_2 are close enough """ - - assert len(array_1) == len(array_2) - for reference_array, test_array in zip(array_1, array_2): - assert len(reference_array) == len(test_array) - assert isclose(reference_array, test_array).all() - class TestPipelinesTeamT54A: """ A class that contains all the unit tests for the PipelineTeamT54A class.""" @@ -77,11 +68,11 @@ def test_outputs(): pipeline = PipelineTeamT54A() # 1 - 1 subject outputs pipeline.subject_list = ['001'] - helpers.test_pipeline_outputs(pipeline, [0, 9*4*1, 5*2*1, 8*2*2 + 4, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 9*4*1, 5*2*1, 8*2*2 + 4*2, 18]) # 2 - 4 subjects outputs pipeline.subject_list = ['001', '002', '003', '004'] - helpers.test_pipeline_outputs(pipeline, [0, 9*4*4, 5*2*4, 8*2*2 + 4, 18]) + helpers.test_pipeline_outputs(pipeline, [0, 9*4*4, 5*2*4, 8*2*2 + 4*2, 18]) @staticmethod @mark.unit_test @@ -101,7 +92,7 @@ def test_subject_information(): bunch = info_missed[0] assert isinstance(bunch, Bunch) assert bunch.conditions == ['trial', 'gain', 'loss', 'difficulty', 'response', 'missed'] - compare_float_2d_arrays(bunch.onsets, [ + helpers.compare_float_2d_arrays(bunch.onsets, [ [4.071, 11.834, 27.535, 36.435], [4.071, 11.834, 27.535, 36.435], [4.071, 11.834, 27.535, 36.435], @@ -109,7 +100,7 @@ def test_subject_information(): [6.459, 14.123, 29.615, 38.723], [19.535] ]) - compare_float_2d_arrays(bunch.durations, [ + helpers.compare_float_2d_arrays(bunch.durations, [ [2.388, 2.289, 2.08, 2.288], [2.388, 2.289, 2.08, 2.288], [2.388, 2.289, 2.08, 2.288], @@ -117,7 +108,7 @@ def test_subject_information(): [0.0, 0.0, 0.0, 0.0], [0.0] ]) - compare_float_2d_arrays(bunch.amplitudes, [ + helpers.compare_float_2d_arrays(bunch.amplitudes, [ [1.0, 1.0, 1.0, 1.0], [14.0, 34.0, 10.0, 16.0], [6.0, 14.0, 15.0, 17.0], @@ -125,35 +116,35 @@ def test_subject_information(): [1.0, 1.0, 1.0, 1.0], [1.0] ]) - assert bunch.regressor_names == None - assert bunch.regressors == None + assert bunch.regressor_names is None + assert bunch.regressors is None bunch = info_ok[0] assert isinstance(bunch, Bunch) assert bunch.conditions == ['trial', 'gain', 'loss', 'difficulty', 'response'] - compare_float_2d_arrays(bunch.onsets, [ + helpers.compare_float_2d_arrays(bunch.onsets, [ [4.071, 11.834, 27.535, 36.435], [4.071, 11.834, 27.535, 36.435], [4.071, 11.834, 27.535, 36.435], [4.071, 11.834, 27.535, 36.435], [6.459, 14.123, 29.615, 38.723] ]) - compare_float_2d_arrays(bunch.durations, [ + helpers.compare_float_2d_arrays(bunch.durations, [ [2.388, 2.289, 2.08, 2.288], [2.388, 2.289, 2.08, 2.288], [2.388, 2.289, 2.08, 2.288], [2.388, 2.289, 2.08, 2.288], [0.0, 0.0, 0.0, 0.0] ]) - compare_float_2d_arrays(bunch.amplitudes, [ + helpers.compare_float_2d_arrays(bunch.amplitudes, [ [1.0, 1.0, 1.0, 1.0], [14.0, 34.0, 10.0, 16.0], [6.0, 14.0, 15.0, 17.0], [1.0, 3.0, 10.0, 9.0], [1.0, 1.0, 1.0, 1.0] ]) - assert bunch.regressor_names == None - assert bunch.regressors == None + assert bunch.regressor_names is None + assert bunch.regressors is None @staticmethod @mark.unit_test @@ -195,10 +186,10 @@ def test_two_sample_t_test_regressors(): ['002', '004'], # equalIndifference group ['001', '002', '003', '004'] # all subjects ) - assert regressors == dict( - equalRange = [1, 0, 1, 0], - equalIndifference = [0, 1, 0, 1] - ) + assert regressors == { + 'equalRange': [1, 0, 1, 0], + 'equalIndifference': [0, 1, 0, 1] + } assert groups == [1, 2, 1, 2] @staticmethod diff --git a/tests/pipelines/test_team_U26C.py b/tests/pipelines/test_team_U26C.py new file mode 100644 index 00000000..c670b520 --- /dev/null +++ b/tests/pipelines/test_team_U26C.py @@ -0,0 +1,142 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Tests of the 'narps_open.pipelines.team_U26C' module. + +Launch this test with PyTest + +Usage: +====== + pytest -q test_team_U26C.py + pytest -q test_team_U26C.py -k +""" +from os import mkdir +from os.path import join, exists +from shutil import rmtree +from filecmp import cmp + +from pytest import helpers, mark, fixture +from numpy import isclose +from nipype import Workflow +from nipype.interfaces.base import Bunch + +from narps_open.utils.configuration import Configuration +from narps_open.pipelines.team_U26C import PipelineTeamU26C + +TEMPORARY_DIR = join(Configuration()['directories']['test_runs'], 'test_U26C') + +@fixture +def remove_test_dir(): + """ A fixture to remove temporary directory created by tests """ + + rmtree(TEMPORARY_DIR, ignore_errors = True) + mkdir(TEMPORARY_DIR) + yield # test runs here + #rmtree(TEMPORARY_DIR, ignore_errors = True) + +def compare_float_2d_arrays(array_1, array_2): + """ Assert array_1 and array_2 are close enough """ + + assert len(array_1) == len(array_2) + for reference_array, test_array in zip(array_1, array_2): + assert len(reference_array) == len(test_array) + assert isclose(reference_array, test_array).all() + +class TestPipelinesTeamU26C: + """ A class that contains all the unit tests for the PipelineTeamU26C class.""" + + @staticmethod + @mark.unit_test + def test_create(): + """ Test the creation of a PipelineTeamU26C object """ + + pipeline = PipelineTeamU26C() + + # 1 - check the parameters + assert pipeline.fwhm == 5.0 + assert pipeline.team_id == 'U26C' + + # 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 PipelineTeamU26C object """ + pipeline = PipelineTeamU26C() + # 1 - 1 subject outputs + pipeline.subject_list = ['001'] + helpers.test_pipeline_outputs(pipeline, [0, 0, 7, 8*3*2 + 5*3, 18]) + + # 2 - 4 subjects outputs + pipeline.subject_list = ['001', '002', '003', '004'] + helpers.test_pipeline_outputs(pipeline, [0, 0, 28, 8*3*2 + 5*3, 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 = PipelineTeamU26C.get_subject_information([test_file, test_file]) + + # Compare bunches to expected + bunch = info[0] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['gamble_run1'] + compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 19.535, 27.535, 36.435]]) + compare_float_2d_arrays(bunch.durations, [[4.0, 4.0, 4.0, 4.0, 4.0]]) + assert bunch.amplitudes == None + assert bunch.tmod == None + assert bunch.pmod[0].name == ['gain_run1', 'loss_run1'] + assert bunch.pmod[0].poly == [1, 1] + 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 == None + assert bunch.regressors == None + + bunch = info[1] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['gamble_run2'] + compare_float_2d_arrays(bunch.onsets, [[4.071, 11.834, 19.535, 27.535, 36.435]]) + compare_float_2d_arrays(bunch.durations, [[4.0, 4.0, 4.0, 4.0, 4.0]]) + assert bunch.amplitudes == None + assert bunch.tmod == None + assert bunch.pmod[0].name == ['gain_run2', 'loss_run2'] + assert bunch.pmod[0].poly == [1, 1] + 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 == None + assert bunch.regressors == None + + @staticmethod + @mark.unit_test + def test_confounds_file(remove_test_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_U26C', 'confounds.tsv') + + # Get new confounds file + PipelineTeamU26C.get_confounds_file(confounds_file, 'sid', 'rid', TEMPORARY_DIR) + + # Check confounds file was created + created_confounds_file = join( + TEMPORARY_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 PipelineTeamU26C and compare results """ + helpers.test_pipeline_evaluation('U26C') diff --git a/tests/pipelines/test_team_X19V.py b/tests/pipelines/test_team_X19V.py new file mode 100644 index 00000000..e0a1a7f3 --- /dev/null +++ b/tests/pipelines/test_team_X19V.py @@ -0,0 +1,138 @@ +#!/usr/bin/python +# coding: utf-8 + +""" Tests of the 'narps_open.pipelines.team_X19V' module. + +Launch this test with PyTest + +Usage: +====== + pytest -q test_team_X19V.py + pytest -q test_team_X19V.py -k +""" +from os import mkdir +from os.path import join, exists +from shutil import rmtree +from filecmp import cmp + +from pytest import helpers, mark, fixture +from numpy import isclose +from nipype import Workflow +from nipype.interfaces.base import Bunch + +from narps_open.utils.configuration import Configuration +from narps_open.pipelines.team_X19V import PipelineTeamX19V + +TEMPORARY_DIR = join(Configuration()['directories']['test_runs'], 'test_X19V') + +@fixture +def remove_test_dir(): + """ A fixture to remove temporary directory created by tests """ + + rmtree(TEMPORARY_DIR, ignore_errors = True) + mkdir(TEMPORARY_DIR) + yield # test runs here + #rmtree(TEMPORARY_DIR, ignore_errors = True) + +def compare_float_2d_arrays(array_1, array_2): + """ Assert array_1 and array_2 are close enough """ + + assert len(array_1) == len(array_2) + for reference_array, test_array in zip(array_1, array_2): + assert len(reference_array) == len(test_array) + assert isclose(reference_array, test_array).all() + +class TestPipelinesTeamX19V: + """ A class that contains all the unit tests for the PipelineTeamX19V class.""" + + @staticmethod + @mark.unit_test + def test_create(): + """ Test the creation of a PipelineTeamX19V object """ + + pipeline = PipelineTeamX19V() + + # 1 - check the parameters + assert pipeline.fwhm == 5.0 + assert pipeline.team_id == 'X19V' + + # 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 PipelineTeamX19V object """ + + pipeline = PipelineTeamX19V() + + # 1 - 1 subject outputs + pipeline.subject_list = ['001'] + helpers.test_pipeline_outputs(pipeline, [0, 4*1 + 4*4*4*1, 4*4*1 + 4*1, 8*4*2 + 4*4, 18]) + + # 2 - 4 subjects outputs + pipeline.subject_list = ['001', '002', '003', '004'] + helpers.test_pipeline_outputs(pipeline, [0, 4*4 + 4*4*4*4, 4*4*4 + 4*4, 8*4*2 + 4*4, 18]) + + @staticmethod + @mark.unit_test + def test_subject_information(): + """ Test the get_subject_information method """ + + # Get test files + test_file = join(Configuration()['directories']['test_data'], 'pipelines', 'events.tsv') + + # Prepare several scenarii + info_ok = PipelineTeamX19V.get_subject_information(test_file) + + # Compare bunches to expected + bunch = info_ok[0] + assert isinstance(bunch, Bunch) + assert bunch.conditions == ['trial', 'gain', 'loss'] + 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]]) + 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]]) + compare_float_2d_arrays(bunch.amplitudes, [ + [1.0, 1.0, 1.0, 1.0, 1.0], + [-8.4, 11.6, 15.6, -12.4, -6.4], + [-8.2, -0.2, 4.8, 0.8, 2.8]]) + assert bunch.regressor_names == None + assert bunch.regressors == None + + @staticmethod + @mark.unit_test + def test_confounds_file(remove_test_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_X19V', 'confounds.tsv') + + # Get new confounds file + PipelineTeamX19V.get_confounds_file(confounds_file, 'sid', 'rid', TEMPORARY_DIR) + + # Check confounds file was created + created_confounds_file = join( + TEMPORARY_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 PipelineTeamX19V and compare results """ + helpers.test_pipeline_evaluation('X19V') diff --git a/tests/test_conftest.py b/tests/test_conftest.py index a9cafbb2..4ea92fb5 100644 --- a/tests/test_conftest.py +++ b/tests/test_conftest.py @@ -25,17 +25,16 @@ from narps_open.utils.configuration import Configuration from narps_open.runner import PipelineRunner from narps_open.pipelines import Pipeline -from narps_open.data.results import ResultsCollection TEST_DIR = abspath(join(Configuration()['directories']['test_runs'], 'test_conftest')) @fixture def set_test_directory(scope = 'function'): + """ A fixture to remove temporary directory created by tests """ + rmtree(TEST_DIR, ignore_errors = True) makedirs(TEST_DIR, exist_ok = True) - yield - # Comment this line for debugging rmtree(TEST_DIR, ignore_errors = True) @@ -234,7 +233,6 @@ def get_file_urls(self): def download(self): """ Download the collection, file by file. """ - pass class TestConftest: """ A class that contains all the unit tests for the conftest module.""" @@ -271,7 +269,7 @@ def test_test_outputs(set_test_directory): # Test pipeline pipeline = MockupPipeline() pipeline.subject_list = ['001', '002'] - + # Wrong length for nb_of_outputs with raises(AssertionError): helpers.test_pipeline_outputs(pipeline, [1,2,3]) @@ -371,21 +369,21 @@ def test_test_pipeline_execution(mocker, set_test_directory): with open(join(TEST_DIR, 'test_conftest.txt'), 'r', encoding = 'utf-8') as file: assert file.readline() == '0\n' # First exec of preprocessing creates an exception (execution counter == 1) - assert file.readline() == f'TestConftest_preprocessing_workflow 4 1\n' + assert file.readline() == 'TestConftest_preprocessing_workflow 4 1\n' # Relaunching the workflow # Preprocessing files won't be created(execution counter == 2) - assert file.readline() == f'TestConftest_preprocessing_workflow 4 2\n' - assert file.readline() == f'TestConftest_run_level_workflow 4 3\n' - assert file.readline() == f'TestConftest_subject_level_workflow 4 4\n' + assert file.readline() == 'TestConftest_preprocessing_workflow 4 2\n' + assert file.readline() == 'TestConftest_run_level_workflow 4 3\n' + assert file.readline() == 'TestConftest_subject_level_workflow 4 4\n' # Relaunching the workflow # Everything's fine - assert file.readline() == f'TestConftest_preprocessing_workflow 4 5\n' - assert file.readline() == f'TestConftest_run_level_workflow 4 6\n' - assert file.readline() == f'TestConftest_subject_level_workflow 4 7\n' - assert file.readline() == f'TestConftest_preprocessing_workflow 3 8\n' - assert file.readline() == f'TestConftest_run_level_workflow 3 9\n' - assert file.readline() == f'TestConftest_subject_level_workflow 3 10\n' - assert file.readline() == f'TestConftest_group_level_workflow 7 11' + assert file.readline() == 'TestConftest_preprocessing_workflow 4 5\n' + assert file.readline() == 'TestConftest_run_level_workflow 4 6\n' + assert file.readline() == 'TestConftest_subject_level_workflow 4 7\n' + assert file.readline() == 'TestConftest_preprocessing_workflow 3 8\n' + assert file.readline() == 'TestConftest_run_level_workflow 3 9\n' + assert file.readline() == 'TestConftest_subject_level_workflow 3 10\n' + assert file.readline() == 'TestConftest_group_level_workflow 7 11' @staticmethod @mark.unit_test diff --git a/tests/test_data/pipelines/team_U26C/confounds.tsv b/tests/test_data/pipelines/team_U26C/confounds.tsv new file mode 100644 index 00000000..05925432 --- /dev/null +++ b/tests/test_data/pipelines/team_U26C/confounds.tsv @@ -0,0 +1,3 @@ +6551.281999999999 6476.4653 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +6484.7285 6473.4890000000005 -0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 -0.00996895 -0.0313444 -3.00931e-06 0.00132687 -0.000384193 -0.00016819 +6441.5337 6485.7256 -2.56954e-05 -0.00923735 0.0549667 0.000997278 -0.00019745 -0.000398988 0.009943254600000001 0.022107050000000003 0.05496970931 -0.00032959199999999986 0.000186743 -0.00023079800000000002 diff --git a/tests/test_data/pipelines/team_X19V/confounds.tsv b/tests/test_data/pipelines/team_X19V/confounds.tsv new file mode 100644 index 00000000..cf63c178 --- /dev/null +++ b/tests/test_data/pipelines/team_X19V/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 diff --git a/tests/test_runner.py b/tests/test_runner.py index bb2a62c3..b72916f9 100644 --- a/tests/test_runner.py +++ b/tests/test_runner.py @@ -142,18 +142,23 @@ def __init__(self): super().__init__() def get_preprocessing(self): + """ Return a preprocessing workflow with wrong type """ return 'Wrong_workflow_type' def get_run_level_analysis(self): + """ Return a run level analysis workflow """ return None def get_subject_level_analysis(self): + """ Return a subject level analysis workflow """ return None def get_group_level_analysis(self): + """ Return a group level analysis workflow """ return None def get_hypotheses_outputs(self): + """ Return hypotheses """ return None class MockupWrongPipeline2(Pipeline): @@ -163,18 +168,23 @@ def __init__(self): super().__init__() def get_preprocessing(self): + """ Return a preprocessing workflow list with wrong types inside """ return ['Wrong_workflow_type', 'Wrong_workflow_type'] def get_run_level_analysis(self): + """ Return a run level analysis workflow """ return None def get_subject_level_analysis(self): + """ Return a subject level analysis workflow """ return None def get_group_level_analysis(self): + """ Return a group level analysis workflow """ return None def get_hypotheses_outputs(self): + """ Return hypotheses """ return None class TestPipelineRunner: