level1_volume.py

#!/usr/bin/env python

# A minimal nipype script to demonstrate first level task fMRI analysis on volume data

import sys
import os
import argparse
import numpy as np
import pandas as pd
import json
import nipype
import nipype.pipeline as pe

import nipype.interfaces.fsl as fsl
import nipype.algorithms.modelgen as modelgen
import nipype.interfaces.io as io
import nipype.interfaces.utility as utility
from nipype.interfaces.base import Bunch
from pathlib import Path

# from nipype import config
# config.enable_debug_mode()
# config.stop_on_first_crash = True

parser = argparse.ArgumentParser(description='First level models for volumetric fmriprep data')
parser.add_argument('bids_dir', help='Full path to the directory with the input dataset '
                    'formatted according to the BIDS standard.')
parser.add_argument('fmriprep_dir', help='Full path to the fmriprep output directory, containing subject level folders')
parser.add_argument('output_dir', help='Full path to the base directory for outputs')
parser.add_argument('tcontrast_file', help='Full path to T contrast specification')
parser.add_argument('participant', help='Participant label to process', type=str)
parser.add_argument('task', help='Task label to process', type=str)
parser.add_argument('run', help='Run label to process', type=str)
parser.add_argument('--session', help='Session label to process',
                   default=None, type=str)
parser.add_argument('--fmriprep_confounds', help='Additional fmriprep confound columns to include.',
                    nargs='*', type=str)
parser.add_argument('--n_cpus','--n-cpus', help='Number of CPUs/cores available to use.',
                   default=1, type=int)
parser.add_argument('--fwhm', help='Smoothing FWHM in mm',
                   default=5.0, type=float)
parser.add_argument('--hpfilter', help='Highpass filter length in seconds',
                   default=120.0, type=float)
parser.add_argument('--drop_trs', help='Number of initial TRs to drop',
                   default=0, type=int)
parser.add_argument('--motion6', help='Include 6 motion regressors as confounds',
                    action='store_true', default=False)
parser.add_argument('--motion6dt', help='Include the first derivative of the 6 motion regressors as confounds',
                    action='store_true', default=False)
parser.add_argument('--space', help='fmriprep output space to use',
                    default='MNI152NLin2009cAsym', type=str)
args = parser.parse_args()

# set up variables used through out the script
subj_id = args.participant
task = args.task
run = args.run
space = args.space # fmriprep output space
# path to processed fmriprep data
fmriprep_path = args.fmriprep_dir
# path to BIDS data
bids_path = args.bids_dir
#path to save results
out_path = args.output_dir

base_dir = f'{out_path}/sub-{subj_id}/'
Path(base_dir).mkdir(parents=True, exist_ok=True)

if args.session is not None:
    session_path = f'ses-{args.session}'
    session_key = f'ses-{args.session}_'
else:
    session_path = ''
    session_key = ''

# get TR
with open(f'{bids_path}/sub-{subj_id}/{session_path}/func/sub-{subj_id}_{session_key}task-{task}_run-{run}_bold.json', 'r') as fp:
    bold_meta = json.load(fp)

# Generate FEAT design files
design = pe.Node(interface = fsl.Level1Design(), name='design', iterfield = [''])
design.inputs.interscan_interval = bold_meta['RepetitionTime']
design.inputs.bases = {'dgamma':{'derivs': False}}
design.inputs.model_serial_correlations = True

# create t contrast list
tcontrasts = pd.read_csv(args.tcontrast_file, sep='\t')
condition_names = list(tcontrasts.columns)[1::]

contrasts = []
for idx, row in tcontrasts.iterrows():
    contrasts.append((row['contrast'], 'T', condition_names, list(row[condition_names])))
design.inputs.contrasts = contrasts


# load confounds
# all columns are demeaned by FSL
confound_names = []
confounds = pd.read_csv(f'{fmriprep_path}/sub-{subj_id}/{session_path}/func/sub-{subj_id}_{session_key}task-{task}_run-{run}_desc-confounds_regressors.tsv',
                            sep='\t')
if args.fmriprep_confounds is not None:
    for concol in args.fmriprep_confounds:
        if concol not in list(confounds.columns):
            raise Exception(f'Confound {concol} not found')
    
    confound_names += args.fmriprep_confounds

confounds = confounds.iloc[args.drop_trs::]
if args.motion6:
    confound_names += ['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z']
if args.motion6dt:
    confound_names += ['trans_x_derivative1', 'trans_y_derivative1', 'trans_z_derivative1', 'rot_x_derivative1', 'rot_y_derivative1', 'rot_z_derivative1']
confound_regressors = [list(confounds[x]) for x in confound_names]

# Make FEAT model
event_file = f'{bids_path}/sub-{subj_id}/{session_path}/func/sub-{subj_id}_{session_key}task-{task}_run-{run}_events.tsv'
if not os.path.exists(event_file):
    print(f'Could not find the run-specific event file {event_file}. Checking for a top-level file.')
    event_file = f'{bids_path}/task-{task}_events.tsv'

if not os.path.exists(event_file):
    raise Exception('Cannot find an _events.tsv file')

trials = pd.read_csv(event_file, sep='\t')

# type checking
assert np.issubdtype(trials['onset'].dtype, np.number), 'Event onsets must be numeric'
assert np.issubdtype(trials['duration'].dtype, np.number), 'Event durations must be numeric'
assert trials.shape[0] == trials[['onset', 'duration', 'trial_type']].dropna().shape[0], 'Event files cannot contain missing values'

trials['onset'] = trials['onset'] - args.drop_trs * bold_meta['RepetitionTime']
trials['trial_type'] = trials['trial_type'].astype(str) # needed in case of numeric trial_type
groups = trials.groupby('trial_type')
onsets = []
durations = []
# check the design
if not set(condition_names).issubset(set(groups.groups.keys())):
    raise Exception(f'Missing conditions in {event_file}')
# keep the correct condition order
for g in condition_names:
    condition_onsets = groups.get_group(g)
    onsets += [list(condition_onsets['onset'])]
    durations += [list(condition_onsets['duration'])]

model_spec = pe.Node(interface = modelgen.SpecifyModel(), name='model_spec', iterfield = [''])
model_spec.inputs.subject_info = [Bunch(conditions=condition_names, onsets=onsets, durations=durations,
                                    regressor_names = confound_names, regressors=confound_regressors)]

model_spec.inputs.input_units = 'secs'
model_spec.inputs.high_pass_filter_cutoff = args.hpfilter
model_spec.inputs.time_repetition = bold_meta['RepetitionTime']

fsl_FEATModel = pe.Node(interface = fsl.FEATModel(), name='fsl_FEATModel', iterfield = [''])

# smooth
fsl_smooth = pe.Node(interface = fsl.SpatialFilter(), name='fsl_smooth', iterfield = [''])
fsl_smooth.inputs.operation = 'mean'
fsl_smooth.inputs.kernel_shape = 'gauss'
# FSL specifies kernels in sd. Convert from FWHM to SD
fsl_smooth.inputs.kernel_size = args.fwhm / (2.0 * (np.sqrt(2.0 * np.log(2.0))))

# drop pre steady state volumes from temporary NIFTI
drop_volumes = pe.Node(interface = fsl.ExtractROI(), name='drop_volumes', iterfield = [''])
drop_volumes.inputs.t_min = args.drop_trs
drop_volumes.inputs.t_size = -1

# brain mask
apply_mask = pe.Node(interface = fsl.ApplyMask(), name='apply_mask', iterfield = [''])

# intensity normalization
fsl_tmean = pe.Node(interface = fsl.ImageStats(op_string='-M'), name='fsl_tmean' )

fn_Scaling = pe.Node(interface = utility.wrappers.Function(input_names=['mean'], output_names=['scaling_op']), name='fn_Scaling')
fn_Scaling.inputs.function_str = "def f(mean): return('-mul ' + str(10000.0/float(mean)))"
fsl_IntNorm = pe.Node(interface = fsl.ImageMaths(), name='fsl_IntNorm')

# Save the prefiltered mean
fsl_MeanImage = pe.Node(interface = fsl.MeanImage(), name='fsl_MeanImage', iterfield = [''])
fsl_MeanImage.inputs.dimension = 'T'

# temporal filtering of the NIFTI
fsl_TemporalFilter = pe.Node(interface = fsl.TemporalFilter(), name='fsl_TemporalFilter', iterfield = [''])
fsl_TemporalFilter.inputs.highpass_sigma = 0.5 * args.hpfilter / bold_meta['RepetitionTime']

# add the mean back to the filtered data
fsl_MultiImageMaths = pe.Node(interface = fsl.MultiImageMaths(), name='fsl_MultiImageMaths', iterfield = [''])
fsl_MultiImageMaths.inputs.op_string = '-add %s'


# GLS model fitting
filmgls_volume = pe.Node(interface = fsl.FILMGLS(), name='filmgls_volume', iterfield = [''])
filmgls_volume.inputs.mode = 'volumetric'
filmgls_volume.inputs.smooth_autocorr = True
filmgls_volume.inputs.mask_size = 5
filmgls_volume.inputs.threshold = 1000
filmgls_volume.inputs.output_type = 'NIFTI_GZ'

#Flexibly collect data from disk to feed into workflows.
select_func = pe.Node(
    io.SelectFiles(
        templates={'out_bold': f'sub-{subj_id}_{session_key}task-{task}_run-{run}_space-{space}_desc-preproc_bold.nii.gz',
                'out_mask': f'sub-{subj_id}_{session_key}task-{task}_run-{run}_space-{space}_desc-brain_mask.nii.gz'}
        ),
    name = 'select_func'
    )
select_func.inputs.base_directory = f'{fmriprep_path}/sub-{subj_id}/{session_path}/func'

# Select the files we want to keep
select_pe_volume = pe.Node(interface = utility.Select(), name='select_pe_volume', iterfield = [''])
select_cope_volume = pe.Node(interface = utility.Select(), name='select_cope_volume', iterfield = [''])
select_varcope_volume = pe.Node(interface = utility.Select(), name='select_varcope_volume', iterfield = [''])
select_zstat_volume = pe.Node(interface = utility.Select(), name='select_zstat_volume', iterfield = [''])

#Basic interface class generates identity mappings
pe_list = pe.Node(utility.IdentityInterface(fields=['pe_index']), name='pe_list')
pe_list.iterables = ('pe_index', np.arange(len(condition_names)))
n_tstats = len([x[0] for x in design.inputs.contrasts if x[1]=='T'])
cope_list = pe.Node(utility.IdentityInterface(fields=['cope_index']), name='cope_list')
cope_list.iterables = ('cope_index', np.arange(n_tstats))

varcope_list = pe.Node(utility.IdentityInterface(fields=['varcope_index']), name='varcope_list')
varcope_list.iterables = ('varcope_index', np.arange(n_tstats))

zstat_list = pe.Node(utility.IdentityInterface(fields=['zstat_index']), name='zstat_list')
zstat_list.iterables = ('zstat_index', np.arange(n_tstats))

#Generic datasink module to store structured outputs
io_DataSink = pe.Node(interface = io.DataSink(), name='io_DataSink', iterfield = [''])
io_DataSink.inputs.base_directory = base_dir
io_DataSink.inputs.container = f'sub-{subj_id}_{session_key}task-{task}_run-{run}_bold'
io_DataSink.inputs.parameterization = False

#Create a workflow to connect all those nodes
analysisflow = nipype.Workflow(f'sub-{subj_id}_{session_key}task-{task}_run-{run}_level1_wf')
#analysisflow = nipype.Workflow(f'stroop_level1_wf')

# drop and mask
analysisflow.connect(select_func, "out_bold", apply_mask, "in_file")
analysisflow.connect(select_func, "out_mask", apply_mask, "mask_file")
analysisflow.connect(apply_mask, "out_file", drop_volumes, "in_file")

#design
analysisflow.connect(drop_volumes, "roi_file", model_spec, "functional_runs")
analysisflow.connect(model_spec, "session_info", design, "session_info")
analysisflow.connect(design, "fsf_files", fsl_FEATModel, "fsf_file")
analysisflow.connect(design, "ev_files", fsl_FEATModel, "ev_files")
analysisflow.connect(fsl_FEATModel, "design_file", filmgls_volume, "design_file")
analysisflow.connect(fsl_FEATModel, "con_file", filmgls_volume, "tcon_file")
analysisflow.connect(fsl_FEATModel, "fcon_file", filmgls_volume, "fcon_file")

# intensity normalization
analysisflow.connect(drop_volumes, "roi_file", fsl_tmean, "in_file")
analysisflow.connect(fsl_tmean, "out_stat", fn_Scaling, "mean")
analysisflow.connect(fn_Scaling, "scaling_op", fsl_IntNorm, "op_string")
analysisflow.connect(drop_volumes, "roi_file", fsl_IntNorm, "in_file")

# temporal filtering
analysisflow.connect(fsl_IntNorm, "out_file", fsl_TemporalFilter, "in_file")
analysisflow.connect(fsl_IntNorm, "out_file", fsl_MeanImage, "in_file")
analysisflow.connect(fsl_TemporalFilter, "out_file", fsl_MultiImageMaths, "in_file")
analysisflow.connect(fsl_MeanImage, "out_file", fsl_MultiImageMaths, "operand_files")

# smooth
analysisflow.connect(fsl_MultiImageMaths, "out_file", fsl_smooth, "in_file")
analysisflow.connect(fsl_smooth, "out_file", filmgls_volume, "in_file")

## result merge and sink

#DF
analysisflow.connect(filmgls_volume, "dof_file", io_DataSink, "stats")

#cope
analysisflow.connect(cope_list, "cope_index", select_cope_volume, "index")
analysisflow.connect(filmgls_volume, "copes", select_cope_volume, "inlist")
analysisflow.connect(select_cope_volume, "out", io_DataSink, "stats.@cope")

#varcopes
analysisflow.connect(varcope_list, "varcope_index", select_varcope_volume, "index")
analysisflow.connect(filmgls_volume, "varcopes", select_varcope_volume, "inlist")
analysisflow.connect(select_varcope_volume, "out", io_DataSink, "stats.@varcope")

# merge outputs into a single CIFTI dense timeseries
# zstats
analysisflow.connect(zstat_list, "zstat_index", select_zstat_volume, "index")
analysisflow.connect(filmgls_volume, "zstats", select_zstat_volume, "inlist")
analysisflow.connect(select_zstat_volume, "out", io_DataSink, "stats.@zstat")

# design debugging
analysisflow.connect(fsl_FEATModel, "design_file", io_DataSink, "design")
analysisflow.connect(fsl_FEATModel, "design_image", io_DataSink, "design.@design_image")
analysisflow.connect(fsl_FEATModel, "design_cov", io_DataSink, "design.@design_cov")
analysisflow.connect(fsl_FEATModel, "con_file", io_DataSink, "design.@con_file")
analysisflow.connect(fsl_FEATModel, "fcon_file", io_DataSink, "design.@fcon_file")

#Run the workflow
plugin = 'MultiProc' #adjust your desired plugin here
plugin_args = {'n_procs': args.n_cpus} #adjust to your number of cores
#analysisflow.write_graph(graph2use='flat', format='png', simple_form=False)
analysisflow.run(plugin=plugin, plugin_args=plugin_args)