07b_GLM-FEAT-SSTMSIT-hp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "071e7555",
   "metadata": {},
   "outputs": [],
   "source": [
    "# from IPython.core.display import display, HTML\n",
    "# display(HTML(\"<style>.container { width:95% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cb3fee6",
   "metadata": {},
   "source": [
    "## FSL FEAT in nipype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b76e75cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "\n",
    "import nipype\n",
    "import nipype.interfaces.io as nio\n",
    "import nipype.interfaces.fsl as fsl\n",
    "import nipype.interfaces.ants as ants\n",
    "import nipype.pipeline.engine as pe\n",
    "import nipype.interfaces.utility as util\n",
    "import nipype.algorithms.modelgen as model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1339730",
   "metadata": {},
   "source": [
    "## Which ses and task?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12968288",
   "metadata": {},
   "outputs": [],
   "source": [
    "# general set-up\n",
    "project_folder = '/home/Public/trondheim'\n",
    "work_dir = os.path.join(project_folder, 'processing', 'nipype_workflow_folders')\n",
    "\n",
    "t_r = 1.38\n",
    "hpcutoff = 10000.   # in seconds\n",
    "\n",
    "ses = 'sstmsit'\n",
    "task = 'msit'\n",
    "model_n = [1]\n",
    "spaces = ['T1w']   # shouldn't touch this but just in case we _do_ want to go back to MNI....\n",
    "smoothing_fwhm = ['1p5'] #, '4p5']\n",
    "\n",
    "template_brain = '/home/Public/trondheim/sourcedata/templates/mni_icbm152_t1_tal_nlin_asym_09c_brain.nii'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1b95590",
   "metadata": {},
   "outputs": [],
   "source": [
    "# which subs can we run?\n",
    "subs_ = [s.split('-')[1] for s in [x.split('/')[3] for x in sorted(glob.glob(f'../derivatives/susan_smoothed_hp_func/sub-*/ses-sstmsit/func/*task-{task}*run-1*fwhm-1p5*'))]]\n",
    "already_run = [s.split('-')[1] for s in [x.split('/')[4] for x in sorted(glob.glob(f'../derivatives/glm_feat_hp_sct/subject_level_model/sub-*/ses-sstmsit/func/fwhm-*/model-{model_n}/*task-{task}*'))]]\n",
    "to_run = [x for x in subs_ if x not in already_run]\n",
    "print(len(to_run))\n",
    "to_run\n",
    "# to_run = [x for x in to_run if '025' not in x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3aa43c48",
   "metadata": {},
   "outputs": [],
   "source": [
    "if task == 'sst':\n",
    "    to_run = ['002','003','004','005','006','007','008','009','010','011','012','013','014','015','016','017','018',\n",
    "              '019','020','021','023','024','025','026','027','029','031','032','033','034','035','037','038','039','040','041','042']\n",
    "    to_remove = ['004', '010', '008', '019', '013', '027']\n",
    "    subject_ids =  [x for x in to_run if x not in to_remove]\n",
    "\n",
    "# no subject 025 for msit (used wrong response buttons)\n",
    "elif task == 'msit':\n",
    "    to_run = ['002','003','004','005','006','007','008','009','010','011','012','013','014','015','016','017','018',\n",
    "              '019','020','021','023','024','025','026','027','029','031','032','033','034','035','037','038','039','040','041','042']\n",
    "    to_remove = ['025']\n",
    "    subject_ids =  [x for x in to_run if x not in to_remove]\n",
    "    \n",
    "# subject_ids = ['040','041']\n",
    "subject_ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a777778",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(subject_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5324945",
   "metadata": {},
   "source": [
    "## First & Second level models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24b80992",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_session_info(subject_id, run, task='sst', ses='sstmsit', space='T1w', shift=-1.38/2, model_n=0,\n",
    "                    root_dir='/home/Public/trondheim', include_physio=True):\n",
    "    # other space: MNI152NLin2009cAsym\n",
    "    import pandas as pd\n",
    "    import numpy as np\n",
    "    import os\n",
    "    from nipype.interfaces.base import Bunch\n",
    "\n",
    "    if task == 'msit':\n",
    "        \n",
    "        ### files ###\n",
    "        sub = subject_id\n",
    "        event_fn = f'{root_dir}/derivatives/event_files/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_events.tsv'\n",
    "#         func_fn = f'{root_dir}/derivatives/fmriprep/fmriprep/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_space-{space}_desc-preproc_bold.nii.gz'\n",
    "        confounds_fn = f'{root_dir}/derivatives/fmriprep/fmriprep/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-confounds_timeseries.tsv'\n",
    "\n",
    "        events = pd.read_csv(event_fn, sep='\\t', index_col=None)\n",
    "        events['duration'] = 0.001  # stick functions\n",
    "        if model_n == 0:\n",
    "            events = events.loc[events.trial_type.isin(['con', 'simon', 'flanker', 'inc'])]\n",
    "        elif model_n == 1:\n",
    "            events = events.loc[events.trial_type.isin(['response_index', 'response_middle', 'response_ring'])]\n",
    "        elif model_n == 2:\n",
    "            events = events.loc[events.trial_type.isin(['correct', 'incorrect'])]\n",
    "            \n",
    "        events = events[['onset', 'trial_type', 'duration']]\n",
    "        \n",
    "        # slice time correction, nb: shift should be a negative number for STC\n",
    "        events['onset'] += shift  \n",
    "\n",
    "        ### confounds ###\n",
    "        confounds = pd.read_csv(confounds_fn, sep='\\t')\n",
    "#        cosine_cols = [x for x in confounds.columns if 'cos' in x]\n",
    "        include_confounds = ['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z', 'dvars', 'framewise_displacement'] #+ cosine_cols\n",
    "        confounds = confounds[include_confounds].fillna(method='bfill')\n",
    "\n",
    "        # get retroicor\n",
    "        if include_physio:\n",
    "            retroicor_fn = f'{root_dir}/derivatives/retroicor/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-retroicor_regressors.tsv'\n",
    "            if not os.path.exists(retroicor_fn):\n",
    "                ## take first 20 aCompCor components\n",
    "                print(\"No retroicor found, including 20 a_comp_cor components\")\n",
    "                a_comp_cor = pd.read_csv(confounds_fn, sep='\\t')[['a_comp_cor_' + str(x).zfill(2) for x in range(20)]]\n",
    "                confounds = pd.concat([confounds, a_comp_cor], axis=1)\n",
    "            else:\n",
    "                retroicor = pd.read_csv(retroicor_fn, sep='\\t', header=None).iloc[:,:20]  ## 20 components in total\n",
    "                retroicor.columns = ['cardiac_' + str(x) for x in range(6)] + ['respiratory_' + str(x) for x in range(8)] + ['respiratoryxcardiac_' + str(x) for x in range(4)] + ['HRV', 'RVT']\n",
    "                confounds = pd.concat([confounds, retroicor], axis=1)\n",
    "\n",
    "        # save the confounds that we actually include in the GLM to csv, no header or index\n",
    "        confounds_fn = f'{root_dir}/derivatives/confounds/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-confounds_timeseries.txt'\n",
    "        os.makedirs(os.path.dirname(confounds_fn), exist_ok=True)\n",
    "        confounds.to_csv(confounds_fn, sep='\\t', header=False, index=False)\n",
    "\n",
    "        ### Make bunch ###\n",
    "        if model_n == 0:\n",
    "            conditions=['con',\n",
    "                        'simon',\n",
    "                        'flanker',\n",
    "                        'inc'\n",
    "                       ]\n",
    "                        \n",
    "        elif model_n == 1:\n",
    "            conditions = ['response_index',\n",
    "                          'response_middle',\n",
    "                          'response_ring'\n",
    "                         ]\n",
    "        \n",
    "        elif model_n == 2:\n",
    "            conditions = ['correct',\n",
    "                          'incorrect']\n",
    "            \n",
    "        onsets = [events.loc[events.trial_type==trial_type,'onset'].tolist() for trial_type in conditions]\n",
    "        durations = [events.loc[events.trial_type==trial_type,'duration'].tolist() for trial_type in conditions]\n",
    "        amplitudes = [[1] * events.loc[events.trial_type==trial_type,'onset'].shape[0] for trial_type in conditions]\n",
    "#             amplitudes = [events.loc[events.trial_type==trial_type,'modulation'].tolist() for trial_type in conditions]\n",
    "\n",
    "        info = Bunch(conditions=conditions,\n",
    "                     onsets=onsets,\n",
    "                     durations=durations,\n",
    "                     amplitudes=amplitudes)\n",
    "    \n",
    "    elif task == 'sst':\n",
    "        \n",
    "        ### files ###\n",
    "        sub = subject_id\n",
    "        event_fn = f'{root_dir}/derivatives/event_files/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_events.tsv'\n",
    "#         func_fn = f'{root_dir}/derivatives/fmriprep/fmriprep/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_space-{space}_desc-preproc_bold.nii.gz'\n",
    "        confounds_fn = f'{root_dir}/derivatives/fmriprep/fmriprep/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-confounds_timeseries.tsv'\n",
    "\n",
    "        events = pd.read_csv(event_fn, sep='\\t', index_col=None)\n",
    "        events['duration'] = 0.001  # stick functions\n",
    "        if model_n == 0:\n",
    "            events = events.loc[events.trial_type.isin(['ss', 'fs', 'go'])]\n",
    "        elif model_n == 1:\n",
    "            events = events.loc[events.trial_type.isin(['response_left', 'response_right'])]\n",
    "        events = events[['onset', 'trial_type', 'duration']]\n",
    "        \n",
    "        # slice time correction, nb: shift should be a negative number for STC\n",
    "        events['onset'] += shift  \n",
    "\n",
    "        ### confounds ###\n",
    "        confounds = pd.read_csv(confounds_fn, sep='\\t')\n",
    "#        cosine_cols = [x for x in confounds.columns if 'cos' in x]\n",
    "        include_confounds = ['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z', 'dvars', 'framewise_displacement'] #+ cosine_cols\n",
    "        confounds = confounds[include_confounds].fillna(method='bfill')\n",
    "\n",
    "        # get retroicor\n",
    "        if include_physio:\n",
    "            retroicor_fn = f'{root_dir}/derivatives/retroicor/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-retroicor_regressors.tsv'\n",
    "            if not os.path.exists(retroicor_fn):\n",
    "                ## take first 20 aCompCor components\n",
    "                print(\"No retroicor found, including 20 a_comp_cor components\")\n",
    "                a_comp_cor = pd.read_csv(confounds_fn, sep='\\t')[['a_comp_cor_' + str(x).zfill(2) for x in range(20)]]\n",
    "                confounds = pd.concat([confounds, a_comp_cor], axis=1)\n",
    "            else:\n",
    "                retroicor = pd.read_csv(retroicor_fn, sep='\\t', header=None).iloc[:,:20]  ## 20 components in total\n",
    "                retroicor.columns = ['cardiac_' + str(x) for x in range(6)] + ['respiratory_' + str(x) for x in range(8)] + ['respiratoryxcardiac_' + str(x) for x in range(4)] + ['HRV', 'RVT']\n",
    "                confounds = pd.concat([confounds, retroicor], axis=1)\n",
    "\n",
    "        # save the confounds that we actually include in the GLM to csv, no header or index\n",
    "        confounds_fn = f'{root_dir}/derivatives/confounds/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-confounds_timeseries.txt'\n",
    "        os.makedirs(os.path.dirname(confounds_fn), exist_ok=True)\n",
    "        confounds.to_csv(confounds_fn, sep='\\t', header=False, index=False)\n",
    "\n",
    "        ### Make bunch ###\n",
    "        if model_n == 0:\n",
    "            conditions=['fs',\n",
    "                        'ss',\n",
    "                        'go']\n",
    "            \n",
    "        elif model_n ==1:\n",
    "            conditions = [\n",
    "                         'response_left',\n",
    "                         'response_right'\n",
    "                         ]\n",
    "\n",
    "        onsets = [events.loc[events.trial_type==trial_type,'onset'].tolist() for trial_type in conditions]\n",
    "        durations = [events.loc[events.trial_type==trial_type,'duration'].tolist() for trial_type in conditions]\n",
    "        amplitudes = [[1] * events.loc[events.trial_type==trial_type,'onset'].shape[0] for trial_type in conditions]\n",
    "#             amplitudes = [events.loc[events.trial_type==trial_type,'modulation'].tolist() for trial_type in conditions]\n",
    "\n",
    "        info = Bunch(conditions=conditions,\n",
    "                     onsets=onsets,\n",
    "                     durations=durations,\n",
    "                     amplitudes=amplitudes)\n",
    "#         print(info)\n",
    "    \n",
    "    # always return info and the confounds_fn\n",
    "    return info, confounds_fn #, contrasts[0]\n",
    "\n",
    "if ses=='sstmsit' and task =='msit' and model_n == [0]:\n",
    "    contrasts = [('con', 'T', ['con'], [1.0]),\n",
    "                 ('simon', 'T', ['simon'], [1.0]),\n",
    "                 ('flanker', 'T', ['flanker'], [1.0]),\n",
    "                 ('inc', 'T', ['inc'], [1.0]),\n",
    "                 ('simon-con', 'T', ['simon', 'con'], [1, -1]),\n",
    "                 ('flanker-con', 'T', ['flanker', 'con'], [1, -1]),\n",
    "                 ('inc-con', 'T', ['inc', 'con'], [1, -1]),\n",
    "                 ('simon-flanker', 'T', ['simon', 'flanker'], [1, -1]),\n",
    "                 ('inc-simon', 'T', ['inc', 'simon'], [1, -1]),\n",
    "                 ('inc-flanker', 'T', ['inc', 'flanker'], [1, -1])\n",
    "                 ]\n",
    "    \n",
    "elif ses=='sstmsit' and task =='msit' and model_n == [1]:\n",
    "    contrasts = [('response_index', 'T', ['response_index'], [1.0]),\n",
    "                 ('response_middle', 'T', ['response_middle'], [1.0]),\n",
    "                 ('response_ring', 'T', ['response_ring'], [1.0]),\n",
    "                 ('index-middle', 'T', ['response_index', 'response_middle'], [1, -1]),\n",
    "                 ('index-ring', 'T', ['response_index','response_ring'], [1, -1]),\n",
    "                 ('middle-ring', 'T', ['response_middle', 'response_ring'], [1, -1])\n",
    "                ]\n",
    "    \n",
    "elif ses=='sstmsit' and task =='msit' and model_n == [2]:\n",
    "    contrasts = [('correct', 'T', ['correct'], [1.0]),\n",
    "                 ('incorrect', 'T', ['incorrect'], [1.0]),\n",
    "                 ('incorrect-correct', 'T', ['incorrect', 'correct'], [1, -1])\n",
    "                ]\n",
    "\n",
    "elif ses=='sstmsit' and task=='sst' and model_n == [0]:\n",
    "    contrasts = [('fs', 'T', ['fs'], [1.0]),\n",
    "                 ('ss', 'T', ['ss'], [1.0]),\n",
    "                 ('go', 'T', ['go'], [1.0]),\n",
    "                 ('fs-go', 'T', ['fs', 'go'], [1, -1]),\n",
    "                 ('fs-ss', 'T', ['fs', 'ss'], [1, -1]),\n",
    "                 ('ss-go', 'T', ['ss', 'go'], [1, -1]),\n",
    "                 ]\n",
    "    \n",
    "elif ses=='sstmsit' and task=='sst' and model_n == [1]:\n",
    "    contrasts = [('response_left', 'T', ['response_left'], [1.0]),\n",
    "                 ('response_right', 'T', ['response_right'], [1.0]),\n",
    "                 ('left-right', 'T', ['response_left','response_right'], [1,-1])\n",
    "                 ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1aa3823",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_runs_per_sub(subject_id):\n",
    "    # here we can add the runs per subject if some subject miss a run or two\n",
    "    # only one run of sst for subs 040 and 041\n",
    "    ### HASH OUT IF RUNNING THE MSIT ###\n",
    "#     if not subject_id in ['040']:\n",
    "#         runs = [1,2]\n",
    "#     else:\n",
    "#         runs = [1]\n",
    "    \n",
    "    ### HASH out if running sst\n",
    "    runs = [1,2]\n",
    "    \n",
    "    return runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb2e7a15",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "workflow = pe.Workflow(name='feat_level12_msit_hp')\n",
    "workflow.base_dir = os.path.join(project_folder, 'processing', 'nipype_workflow_folders')\n",
    "workflow.config = {\"execution\": {\"crashdump_dir\":os.path.join(project_folder, 'processing', 'crashdumps')}}\n",
    "\n",
    "# identity\n",
    "identity = pe.Node(util.IdentityInterface(fields=['subject_id', 'space', 'model_n', 'smoothing_fwhm']), name='identity')\n",
    "identity.iterables = [('subject_id', subject_ids),\n",
    "                      ('space', spaces),\n",
    "                      ('model_n', model_n),\n",
    "                      ('smoothing_fwhm', smoothing_fwhm)]\n",
    "\n",
    "# selector\n",
    "# if running only one run, copy the mask from run 1 to run 2, 'mask' variable below needs a list, it will crash if only one is found\n",
    "#(terrible way to get around this.. fix at some point.. lazy..)\n",
    "templates = {'smoothed_funcs': os.path.join(project_folder, 'derivatives', 'susan_smoothed_hp_func', \n",
    "                                            'sub-{subject_id}',  f'ses-{ses}', 'func',  \n",
    "                                            f'sub-{{subject_id}}_ses-{ses}_task-{task}_run-*_space-{{space}}*_fwhm-{{smoothing_fwhm}}_tempfilt.nii.gz'),\n",
    "             'mask': os.path.join(project_folder, 'derivatives', 'fmriprep', 'fmriprep', \n",
    "                                  'sub-{subject_id}', f'ses-{ses}', 'func', \n",
    "                                  f'sub-{{subject_id}}_ses-{ses}_task-{task}_run-*_space-{{space}}_desc-brain_mask.nii.gz'),\n",
    "             'composite':  os.path.join(project_folder, 'derivatives', 'fmriprep', 'fmriprep', \n",
    "                                        'sub-{subject_id}', 'anat', \n",
    "                                        'sub-{subject_id}_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5'),\n",
    "             'xfm': os.path.join(project_folder, 'derivatives', 'fmriprep', 'fmriprep', \n",
    "                                  'sub-{subject_id}', f'ses-{ses}', 'func', \n",
    "                                  f'sub-{{subject_id}}_ses-{ses}_task-{task}_run-*_from-scanner_to-T1w_mode-image_xfm.txt'),}\n",
    "selector = pe.Node(nio.SelectFiles(templates), name='selector')\n",
    "\n",
    "workflow.connect(identity, 'subject_id', selector, 'subject_id')\n",
    "workflow.connect(identity, 'space', selector, 'space')\n",
    "workflow.connect(identity, 'smoothing_fwhm', selector, 'smoothing_fwhm')\n",
    "\n",
    "# get run info per sub\n",
    "run_info_getter = pe.Node(util.Function(function=get_runs_per_sub,\n",
    "                                        input_names=['subject_id'],\n",
    "                                        output_names=['runs']), name='run_info_getter')\n",
    "workflow.connect(identity, 'subject_id', run_info_getter, 'subject_id')\n",
    "\n",
    "# session info getter\n",
    "session_info_getter = pe.MapNode(util.Function(function=get_session_info,\n",
    "                                     input_names=['subject_id', 'run', 'task', 'ses', 'space', 'shift', 'model_n'],\n",
    "                                     output_names=['session_info', 'confounds']),\n",
    "                                 iterfield=['run'],\n",
    "                                 name='session_info_getter')\n",
    "session_info_getter.inputs.task=task\n",
    "session_info_getter.inputs.ses=ses\n",
    "session_info_getter.inputs.space='T1w'\n",
    "session_info_getter.inputs.shift = -t_r/2\n",
    "\n",
    "workflow.connect(identity, 'subject_id', session_info_getter, 'subject_id')\n",
    "workflow.connect(run_info_getter, 'runs', session_info_getter, 'run')\n",
    "workflow.connect(identity, 'model_n', session_info_getter, 'model_n')\n",
    "\n",
    "\n",
    "# model setup\n",
    "specifymodel = pe.Node(model.SpecifyModel(), name='specifymodel1')\n",
    "specifymodel.inputs.input_units = 'secs'\n",
    "specifymodel.inputs.time_repetition = t_r\n",
    "specifymodel.inputs.high_pass_filter_cutoff = hpcutoff\n",
    "\n",
    "workflow.connect(session_info_getter, 'session_info', specifymodel, 'subject_info')\n",
    "## old flow (with cosines): immediately connect to specifymodel\n",
    "workflow.connect(selector, 'smoothed_funcs', specifymodel, 'functional_runs')\n",
    "\n",
    "# Level 1 design\n",
    "level1design = pe.Node(interface=fsl.Level1Design(), name=\"level1design\")\n",
    "level1design.inputs.interscan_interval = t_r\n",
    "level1design.inputs.bases = {'dgamma': {'derivs': True}}\n",
    "level1design.inputs.contrasts = contrasts\n",
    "level1design.inputs.model_serial_correlations = True\n",
    "\n",
    "workflow.connect(specifymodel, 'session_info', level1design, 'session_info')\n",
    "# workflow.connect(session_info_getter, 'contrasts', level1design, 'contrasts')\n",
    "\n",
    "# FEAT model\n",
    "modelgen = pe.MapNode(interface=fsl.FEATModel(), iterfield=['ev_files', 'fsf_file', 'args'], name='modelgen')\n",
    "\n",
    "workflow.connect(level1design, 'ev_files', modelgen, 'ev_files')\n",
    "workflow.connect(level1design, 'fsf_files', modelgen, 'fsf_file')\n",
    "workflow.connect(session_info_getter, 'confounds', modelgen, 'args')   # add confounds here\n",
    "\n",
    "\n",
    "# FILM GLS\n",
    "iterfield = ['design_file', 'in_file', 'tcon_file']\n",
    "modelestimate = pe.MapNode(interface=fsl.FILMGLS(smooth_autocorr=True,\n",
    "                                                 mask_size = 5,\n",
    "                                                 threshold = 1000),  # \n",
    "                          name='modelestimate',\n",
    "                          iterfield=iterfield,\n",
    "                          mem_gb=10)\n",
    "\n",
    "# ## old flow (with cosines): immediately connect to modelestimate\n",
    "workflow.connect(selector, 'smoothed_funcs', modelestimate, 'in_file')\n",
    "# # new flow: connect high-passed data\n",
    "# workflow.connect(highpass, 'out_file', modelestimate, 'in_file')\n",
    "\n",
    "\n",
    "workflow.connect(modelgen, 'design_file', modelestimate, 'design_file')\n",
    "workflow.connect(modelgen, 'con_file', modelestimate, 'tcon_file')\n",
    "\n",
    "\n",
    "#### Fixed effects\n",
    "# merge copes, varcopes\n",
    "copemerge = pe.MapNode(\n",
    "    interface=fsl.Merge(dimension='t'),\n",
    "    iterfield=['in_files'],\n",
    "    name=\"copemerge\")\n",
    "\n",
    "varcopemerge = pe.MapNode(\n",
    "    interface=fsl.Merge(dimension='t'),\n",
    "    iterfield=['in_files'],\n",
    "    name=\"varcopemerge\")\n",
    "\n",
    "maskemerge = pe.MapNode(interface=fsl.Merge(dimension='t'),\n",
    "                       iterfield=['in_files'],\n",
    "                       name=\"maskemerge\")\n",
    "\n",
    "def sort_copes(files):\n",
    "    numelements = len(files[0])\n",
    "    outfiles = []\n",
    "    for i in range(numelements):\n",
    "        outfiles.insert(i,[])\n",
    "        for j, elements in enumerate(files):\n",
    "            outfiles[i].append(elements[i])\n",
    "    return outfiles\n",
    "\n",
    "workflow.connect(modelestimate, ('copes',sort_copes), copemerge, 'in_files')\n",
    "workflow.connect(modelestimate, ('varcopes',sort_copes), varcopemerge, 'in_files')\n",
    "\n",
    "\n",
    "level2model = pe.Node(interface=fsl.L2Model(), name='l2model')\n",
    "def num_copes(files):\n",
    "    return len(files)\n",
    "workflow.connect(modelestimate, ('copes',num_copes), level2model, 'num_copes')\n",
    "\n",
    "\n",
    "pickfirst = lambda x: x[0]\n",
    "flameo = pe.MapNode(\n",
    "    interface=fsl.FLAMEO(run_mode='fe'),\n",
    "    name=\"flameo\",\n",
    "    iterfield=['cope_file', 'var_cope_file'])\n",
    "\n",
    "workflow.connect([\n",
    "    (selector, flameo, [(('mask', pickfirst), 'mask_file')]),\n",
    "    (copemerge, flameo, [('merged_file', 'cope_file')]),\n",
    "    (varcopemerge, flameo, [('merged_file', 'var_cope_file')]),\n",
    "    (level2model, flameo, [('design_mat', 'design_file'),\n",
    "                           ('design_con', 't_con_file'), \n",
    "                           ('design_grp', 'cov_split_file')]),\n",
    "])\n",
    "\n",
    "## Warp to MNI-space\n",
    "def warp_files(copes, varcopes, zstats, tdofs, mat, template_brain):\n",
    "    # stolen from https://dpaniukov.github.io/2016/07/14/three-level-analysis-with-fsl-and-ants-2.html\n",
    "    import nipype.interfaces.ants as ants\n",
    "\n",
    "    out_copes = []\n",
    "    out_varcopes = []\n",
    "    out_zstats = []\n",
    "    out_tdofs = []\n",
    "    \n",
    "    warp = ants.ApplyTransforms()\n",
    "    warp.inputs.input_image_type = 0\n",
    "    warp.inputs.interpolation = 'Linear'\n",
    "    warp.inputs.invert_transform_flags = [False] #,False]\n",
    "    warp.inputs.reference_image = template_brain\n",
    "    warp.inputs.transforms = mat\n",
    "\n",
    "    if not isinstance(copes, list):\n",
    "        copes = [copes]\n",
    "        varcopes = [varcopes]\n",
    "        zstats = [zstats]\n",
    "        tdofs = [tdofs]\n",
    "    \n",
    "    for cope in copes:\n",
    "        warp.inputs.input_image = cope\n",
    "        res=warp.run()\n",
    "        out_copes.append(str(res.outputs.output_image))\n",
    "\n",
    "    for varcope in varcopes:\n",
    "        warp.inputs.input_image = varcope\n",
    "        res=warp.run()\n",
    "        out_varcopes.append(str(res.outputs.output_image))\n",
    "        \n",
    "    for zstat in zstats:\n",
    "        warp.inputs.input_image = zstat\n",
    "        res=warp.run()\n",
    "        out_zstats.append(str(res.outputs.output_image))\n",
    "    \n",
    "    for tdof in tdofs:\n",
    "        warp.inputs.input_image = tdof\n",
    "        res=warp.run()\n",
    "        out_tdofs.append(str(res.outputs.output_image))\n",
    "\n",
    "    return out_copes, out_varcopes, out_zstats, out_tdofs\n",
    "\n",
    "warpfunc = pe.MapNode(util.Function(input_names=['copes', 'varcopes', 'zstats', 'tdofs', 'mat','template_brain'],\n",
    "                               output_names=['out_copes', 'out_varcopes', 'out_zstats', 'out_tdofs'],\n",
    "                               function=warp_files),\n",
    "                               iterfield=['copes', 'varcopes', 'zstats', 'tdofs'],\n",
    "                  name='warpfunc')\n",
    "\n",
    "warpfunc.inputs.template_brain = template_brain\n",
    "workflow.connect(flameo, 'copes', warpfunc, 'copes')\n",
    "workflow.connect(flameo, 'var_copes', warpfunc, 'varcopes')\n",
    "workflow.connect(flameo, 'zstats', warpfunc, 'zstats')\n",
    "workflow.connect(flameo, 'tdof', warpfunc, 'tdofs')\n",
    "workflow.connect(selector, 'composite', warpfunc, 'mat')\n",
    "\n",
    "#### post-FE analyis\n",
    "# Smoothness estimation\n",
    "smoothestimate = pe.MapNode(fsl.SmoothEstimate(), iterfield=['zstat_file'], name='smoothestimate')\n",
    "\n",
    "workflow.connect(selector, ('mask', pickfirst), smoothestimate, 'mask_file')\n",
    "workflow.connect(flameo, 'zstats', smoothestimate, 'zstat_file')\n",
    "\n",
    "#\n",
    "get_volume = pe.Node(fsl.ImageStats(op_string = '-V'), name='get_volume')\n",
    "workflow.connect(selector, ('mask', pickfirst), get_volume, 'in_file')\n",
    "\n",
    "# Cluster threshold\n",
    "grf_cluster = pe.MapNode(fsl.Cluster(), iterfield=['dlh', 'in_file'], name='grf_cluster')\n",
    "grf_cluster.iterables = [(\"threshold\", [2.3, 3.1])] \n",
    "workflow.connect(smoothestimate, 'dlh', grf_cluster, 'dlh')\n",
    "workflow.connect(flameo, 'zstats', grf_cluster, 'in_file')\n",
    "\n",
    "def convert_volume(input):\n",
    "    return int(input[1])\n",
    "\n",
    "workflow.connect(get_volume, ('out_stat', convert_volume), grf_cluster, 'volume')\n",
    "grf_cluster.inputs.out_threshold_file = True\n",
    "\n",
    "\n",
    "## datasink\n",
    "ds = pe.Node(nio.DataSink(), name='datasink')\n",
    "ds.inputs.base_directory = os.path.join(project_folder, 'derivatives', 'glm_feat_hp_sct')\n",
    "\n",
    "#_model_n_0_smoothing_fwhm_1p5_space_T1w_subject_id_002\n",
    "substitutions = [(f't1w/level2_{stat_type}s/_model_n_{model_n_}_smoothing_fwhm_{fwhm}_space_T1w_subject_id_{sub}/_flameo{contrast_n}/{stat_type}1.nii.gz',\n",
    "                  f'sub-{sub}/ses-{ses}/func/fwhm-{fwhm}/model-{model_n_}/sub-{sub}_ses-{ses}_task-{task}_space-T1w_model-{model_n_}_contrast-{contrast_n}_desc-{stat_type}.nii.gz')\n",
    "                  for sub in subject_ids\n",
    "                  for contrast_n in range(len(contrasts))\n",
    "                  for model_n_ in model_n\n",
    "                  for fwhm in smoothing_fwhm\n",
    "                  for stat_type in ['cope', 'zstat', 'varcope', 'tdof_t']\n",
    "                  ]\n",
    "\n",
    "\n",
    "substitutions += [(f'mni/level2_{stat_type}s/_model_n_{model_n_}_smoothing_fwhm_{fwhm}_space_T1w_subject_id_{sub}/_warpfunc{contrast_n}/{stat_type}1_trans.nii.gz',\n",
    "                   f'sub-{sub}/ses-{ses}/func/fwhm-{fwhm}/model-{model_n_}/sub-{sub}_ses-{ses}_task-{task}_space-MNI152NLin2009cAsym_model-{model_n_}_contrast-{contrast_n}_desc-{stat_type}.nii.gz')\n",
    "                  for sub in subject_ids\n",
    "                  for contrast_n in range(len(contrasts))\n",
    "                  for model_n_ in model_n\n",
    "                  for fwhm in smoothing_fwhm\n",
    "                  for stat_type in ['cope', 'zstat', 'varcope', 'tdof_t']\n",
    "                  ]\n",
    "\n",
    "ds.inputs.substitutions = substitutions\n",
    "\n",
    "\n",
    "## cluster thresholds: leave for now, this is more of a third-level model thing anyway\n",
    "# workflow.connect(grf_cluster, 'threshold_file', ds, 'grf_thresholded_zstats_file')\n",
    "# workflow.connect(grf_cluster, 'locnalmax_txt_file', ds, 'grf_localmax_txt_file')\n",
    "\n",
    "## T1w-space\n",
    "workflow.connect(flameo, 'zstats', ds, 'subject_level_model.t1w.level2_zstats')\n",
    "workflow.connect(flameo, 'copes', ds, 'subject_level_model.t1w.level2_copes')\n",
    "workflow.connect(flameo, 'var_copes', ds, 'subject_level_model.t1w.level2_varcopes')\n",
    "workflow.connect(flameo, 'tdof', ds, 'subject_level_model.t1w.level2_tdof_ts')\n",
    "\n",
    "## MNI-space\n",
    "workflow.connect(warpfunc, 'out_zstats', ds, 'subject_level_model.mni.level2_zstats')\n",
    "workflow.connect(warpfunc, 'out_copes', ds, 'subject_level_model.mni.level2_copes')\n",
    "workflow.connect(warpfunc, 'out_varcopes', ds, 'subject_level_model.mni.level2_varcopes')\n",
    "workflow.connect(warpfunc, 'out_tdofs', ds, 'subject_level_model.mni.level2_tdof_ts')\n",
    "# workflow.connect(flameo, 'tdof', ds, 'level2_tdof')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2ca2d97",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "workflow.run(plugin='MultiProc', plugin_args={'n_procs':12, 'memory_gb': 100})"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py38",
   "language": "python",
   "name": "py38"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}