Usage
For all steps, you can specify what data to process using csv-files. See NHP-BIDS/csv/SubSesRun.csv for an example.
Most steps can in principle also be run on the SurfSara LISA cluster. If you want to do that, you should create a job-file to work with the batch scheduler (see documentation). The precise names and subfolder location of job-files are optional, but it should be a bash script. A template can be found in NHP-BIDS/code/lisa/template_SLURM.sh.
For all steps below, we specify how you would run each processing step on the cluster, but a single job-file can also contain all these steps together in a sequence so that everything runs serially (make sure you reserve enough computing time in your jobs-script!)
NB! To comply with institutional data storage policies, we have made some changes to the organisation of the data. These will now be separated by projects in project-based folders. Inside these folders we maintain the old BIDS conventions. As a consequence, the code for the pipelines has been altered to work with the new path definitions in the code branch nin-dataprotocol. This also affects the LISA cluster and the project disk there. Importantly, it has not been tested yet! For now, Chris Klink also maintains a backup of the old data organisation. The switch will be made completely at some point, but it's a work-in-progress.
Make sure you have the data in a Data_raw folder. The precise method to do this tends to change a bit because it adapts to the most current Spinoza Center data-server solution. At the NIN, Chris Klink (c.klink@nin.knaw.nl) will most likely take care of this export of the data from the scanner and stimulus & eye-tracker computers, but right now (20200201) the workflow is:
- create a
Data_raw/<SUBJECT>/<YYYYMMDD>folder and make subfoldersMRI,Behavior,andEye. - copy the eye-log and behavioral-logs to their respective folders.
- copy the
.nii.gzand corresponding.jsonfiles to theMRIfolder. At the moment the dicom to nifti conversion is done on the SC Flywheel data server and you can download all files directly. Alternatively, you can download just the dicom files (e.g., from the SC ftp server) and do the dicom to nifti conversion yourself. In that case, use thedcm2niixtool and also export the json file:dcm2nnix -o <outputfolder> -b y <dicomfolder>.
Create a copy-to-bids.sh script in the Data_raw/<SUBJECT>/<YYYYMMDD> folder, and run it.
- Base this script on an already existing script from a previous run or simply use the template
NHP-BIDS\code\copy-to-bids_template.sh. - If you have trouble matching imaging data to log-files (e.g., because something went wrong with the filenames) you can either check your written notes, or (more reliably) you need to match up the AcquisitionTime (in the json files, created by dcm2niix) with the behavioral timestamps. If there are no json-files available you can (re)create them with
dcm2niix -b y -o <outputfolder> <location of dicom images>. In addition, the json in the behavioral log folder can tell you under what run number the behavioral session was logged. If you did everything correctly, this should match the run number in the scan's filename. If none of this works, you're out of luck and should log things better next time. -
NB! This is where the new data organisation already comes into play. The new
copy-to-bids.shfile in thenin-dataprotocolbranch of the data works with project folders.
Create or modify a csv file that lists the subject, session, run, and datatype to process (see SubSesRun.csv for an example of how to format this). You can also specify a reference subject (refsubject), but this won't be used until the pre-processing step. The csv files can be kept in the csv directory.
NB! In all workflows except for bids_minimal_processing, the run column specifies which functional runs will be processed. In bids_minimal_processing, all image files are processed by default and the run column specifies which eventlog csv files will be processed. In most cases, you don't need to worry about this since as long as the log-csv for a run is properly created (something to keep in mind when programming new experiments with Tracker) the definition will be the same and csv files will need to be processed for all existing runs anyway.
This workflow will reorient the animal data from sphinx position to regular position so that the axes are correct and show up nicely in viewers like fsleyes. Run ./code/bids_minimal_preprocessing.py from your NHP-BIDS root directory (this file also has instructions in the file header).
- The
--projflag is mandatory and indicates to which project these files belong. There should be a project folder inNHP-BIDS/projectswith the same label. This only hold for thenin-dataprotocolbranch for now but will become the default after testing. - The
--csvflag is mandatory and should point to the csv-file that is formatted as explained above. - The
--ignore_eventsflag causes the workflow to skip processing the eventlog csv files for all runs. If this flag is not used, the event-logs will be processed for all runs specified in the csv file. - example:
clear && ./code/bids_minimal_processing.py --proj <ProjectName> --csv ./code/csv/<SubSesRun.csv> |& tee ./logs/log-minproc.txt - help:
./code/bids_minimal_processing.py --help
- Log on to
lisa.surfsara.nland go to yourNHP-BIDSdirectory. Runsbatch ./code/minproc/<minproc_script>.sh, where <minproc_script.sh> is a correctly formatted job-file. A command like the above should be part of the job-file.
Run ./code/bids_resample_isotropic_workflow.py to resample all volumes to 1.0 mm isotropic voxels
- The
--projflag is mandatory and indicates to which project these files belong. There should be a project folder inNHP-BIDS/projectswith the same label. This only hold for thenin-dataprotocolbranch for now but will become the default after testing. - The
--csvflag is mandatory and should point to the csv-file that is formatted as explained above. - example:
clear && ./code/bids_resample_isotropic_workflow.py --proj <ProjectName> --csv ./code/csv/<SubSesRun.csv> |& tee ./logs/log-resample.txt
- Log on to
lisa.surfsara.nland go to yourNHP-BIDSdirectory. Runsbatch ./code/isoresample/<isoresample_script>.sh, job-file. A command like the above should be part of the job-file.
Run ./code/bids_resample_hires_isotropic_workflow.py if you also want the high-resolution 0.6 mm isotropic anatomical images
- The
--projflag is mandatory and indicates to which project these files belong. There should be a project folder inNHP-BIDS/projectswith the same label. This only hold for thenin-dataprotocolbranch for now but will become the default after testing. - The
--csvflag is mandatory and should point to the csv-file that is formatted as explained above. - example:
clear && ./code/bids_resample_hires_isotropic_workflow.py --proj <ProjectName> --csv ./code/csv/<SubSesRun.csv> |& tee ./logs/log-resample_hiresanat.txt
- Log on to
lisa.surfsara.nl, go to yourNHP-BIDSdirectory, and runsbatch ./code/isoresample/<isoresample_hires_script.sh>.
This workflow ./code/bids_undistort_workflow.py takes the 'topup' scans from the fmap directory for each functional run and undistorts the epi's by calculating the halfway non-linear warp between subsequent scans with opposite phase-encoding direction. The appropriate warp is then applied to the functional run to facilitate alignment.
- This workflow requires the opposite phase-encoding files to exist and will fail if they don't exist.
- The
--projflag is mandatory and indicates to which project these files belong. There should be a project folder inNHP-BIDS/projectswith the same label. This only hold for thenin-dataprotocolbranch for now but will become the default after testing. - The
--csvflag is mandatory and should point to the csv-file that is formatted as explained above. - example:
clear && ./code/bids_undistort_workflow.py --proj <ProjectName> --csv ./code/csv/<SubSesRun.csv> |& tee ./logs/log-undistort.txt
- Log on to
lisa.surfsara.nl, go to yourNHP-BIDSdirectory, and runsbatch ./code/undistort/<undistort_script.sh>.
This workflow ./code/bids_preprocessing_workflow.py does non-linear slice-by-slice alignment, alignment to a functional reference, performs motion correction, and some filtering and smoothing.
- Data is nonlinearly registered to reference volumes that are located in
NHP-BIDS/reference-vols/sub-<subject>/ - This workflow requires references (& warps) for each subject, saved in subject-specific folders in
reference-vols. If you work with a new subject, you should create such a folder. In that case you should follow the filename pattern of the existing references and only replace the subject name. The reference-vols files are too large to maintain in the repo, but they can be found on theVS03-2server and on Dropbox, including Jupyter notebooks with code that can guide you through the creation process.
- The
--projflag is mandatory and indicates to which project these files belong. There should be a project folder inNHP-BIDS/projectswith the same label. This only hold for thenin-dataprotocolbranch for now but will become the default after testing. - The
--csvflag is mandatory and should point to the csv-file that is formatted as explained above. - In this workflow you can add an
refsubjectcolumn to the csv file. The workflow will interpret this column as identifier for which reference folder to use (default is the one named sub-). When the column is omitted, the workflow looks for reference files in the default folder. This can be useful if, for some reason, an animal should be registered against another reference later on (e.g., after reimplantation of a head-post). When used as a header,refsubjectwill have to be specified for each row of the csv file (but you can simply repeat thesubjectentry if you want to use the default for some runs).
NB! Therefsubjectargument specifies which reference-vols folder to get references from. The actual files in this folder should still belong to correct subjects and be named appropriately. Example, you can preprocess files for subject A withreference-vols/sub-A/func/sub-A_ref_etc.nii.gzandreference-vols/sub-A_alternative/func/sub-A_ref_etc.nii.gz, but not withreference-vols/sub-A_alternative/func/sub-B_ref_etc.nii.gz. - The
--undistortflag indicates whether to use the undistorted functional runs (default=True) or the ones without undistortion. - The
--fwhmand--HighPassflags are optional to specify spatial smoothing (mm) and a high-pass filter (s). If not specified, they will be 2.0 mm and 50 s respectively. - example:
clear && ./code/bids_preprocessing_workflow.py --proj <ProjectName> --csv ./code/csv/<SubSesRun.csv> |& tee ./logs/log-preproc.txt
- Log on to
lisa.surfsara.nl, go to yourNHP-BIDSdirectory and runsbatch ./code/lisa/preproc/<preprocess_script.sh>.
NB! It is possible (though rare) that the workflow crashes with a messsage that starts withRuntimeError: Command: convert_xfm -omat ....This is an FSL bug in which a flirt operation creates a hexadecimal matrix file instead of a decimal one. You can fix this with the script./helper_scripts/hex2dec.shand re-run the workflow. - There's an alternative way of doing this on the cluster. If you run things as above, the slice-by-slice alignment will be parallel but the different runs are processed serially. If you want to also do the runs in parallel (making things a lot faster, but requiring more nodes on the cluster) you can use the script
bids_preprocessing_parallel_runs.py. Instructions on how to use it can be found under Parallel preprocessing on LISA. Due to the significant speed increase, this is the preferred method.
- There's now also a way of running this as a parallel process on a local machine (not based on SLURM jobs). Use the
bids_preproc_parallel.pywith an additional--maxprocargument to specify the maximum number of parallel processes. When choosing this, keep the hardware limitations of the current machine in mind and do not request more than you have available. NB! THis approach is currently still untested (20231025). - example:
clear && ./code/bids_preproc_parallel.py --proj <ProjectName> --csv ./code/csv/<SubSesRun.csv> --maxproc 8 |& tee ./logs/log-preproc-par.txt
The motion outlier detection that is part of bids_preprocessing_workflow.py is based on nipype's rapidart algorithm. In our experience, this is sometimes not strict enough. This workflow bids_fslmotionoutliers_workflow.py uses an additional step of FSL's motion outlier detection on the motion corrected time-series and creates merged outlier files that can be used as alternative stricter outlier files when performing statistical analysis. Of course, this means that preprocessing will already have to have been done on any files selected for this workflow. It will output additional files to the ./derivatives/featpreproc/motion_outliers folder containing the FSL confound matrices and a text file of merged identified confound volume indeces in AFNI style. This allows easy usage in later processing phases.
- The
--projflag is mandatory and indicates to which project these files belong. There should be a project folder inNHP-BIDS/projectswith the same label. This only hold for thenin-dataprotocolbranch for now but will become the default after testing. - The
--csvflag is mandatory and should point to the csv-file that is formatted as explained above. - The
--undistortflag indicates whether to use the undistorted functional runs (default=True) or the ones without undistortion.
This workflow (./code/bids_warp2nmt_workflow.py) warps the preprocessed highpassed_files to the common NMT space to facilitate statistical inferences in a common space. The relevant warps are pre-calculated and stored in the reference-vols folder.
- The
--projflag is mandatory and indicates to which project these files belong. There should be a project folder inNHP-BIDS/projectswith the same label. This only hold for thenin-dataprotocolbranch for now but will become the default after testing. - The
--csvflag is mandatory and should point to the csv-file that is formatted as explained above. - The
--undistortflag indicates whether to use the undistorted functional runs (default=True) or the ones without undistortion. - Here you can again specify a refsubject column in the csv.
- example:
clear && ./code/bids_warp2nmt_workflow.py --proj <ProjectName> --csv ./code/csv/warp2nmt/<SubSesRun.csv> |& tee ./logs/log-warp2nmt.txt
- Log on to
lisa.surfsara.nl, go toNHP-BIDSdirectory and runsbatch ./code/lisa/warp2nmt/<SCRIPTNAME>.sh. A command like the above should be part of the job-file. - You can include this step in the fully parallel preprocessing (recommended) with the script
bids_preprocessing_parallel_runs.py. Instructions on how to use it can be found in the section Parallel preprocessing on LISA.
The level1 analysis in bids_modelfit_level1_workflow.py performs a GLM analysis on the preprocessed data of individual runs and outputs statistical results to ./derivatives/modelfit/<contrast>/<registration_space>/level1/mo-<motion_outlier_type>.
NB! All files in the model-fit are expected to be motion corrected and registered to the same reference space!! One way to do get everything to register motion corrected files to a standard space first using the bids_warp2nmt_workflow.py workflow.
- The
--projflag is mandatory and indicates to which project these files belong. There should be a project folder inNHP-BIDS/projectswith the same label. This only hold for thenin-dataprotocolbranch for now but will become the default after testing. - The
--csvflag is mandatory and should point to the csv-file that is formatted as explained above. Here you can again specify a refsubject column in the csv. Make sure this is done correctly as it determines which brainmask file will be used in the analysis. - The
--undistortflag indicates whether to use the undistorted functional runs (default=True) or the ones without undistortion. - The
--contrastsflag is mandatory and depends on the experiment. Incode/contrasts/there are python modules for each set of contrasts. If you want to use the contrasts defined inctcheckerboard.py, for example, passctcheckerboardas the value for the--contrastsparameter. If you create your own contrasts, you just need your function to define the variable namedcontrasts. Since the code assumes that this will be a python module name, you cannot use dashes or spaces in the name. - The optional argument
--RegSpaceallows us to specify in which space the GLM is performed. The default isnmt, which means that the workflow will use the highpassed files that were aligned to theNMTtemplate space as input to the GLM. This puts all subjects in a common reference space, but it only works ifbids_warp2nmt_workflow.pyhas been executed already. The alternative optionnativewill use the preprocessed functional timeseries in their functional space. Be careful with this as there is no explicit check for mismatching functional spaces with multiple entries and statistics will be meaningless if not calculated in a common space for all functional runs. - The optional argument
--resultfldallows you to define the name of the working directory for this analysis. Re-running something with the same name uses existing intermediate results from working directories. Output directories are determined based on the contrast and registration space that are chosen allowing straightforward selection of results for level2 analysis. If you do not include a `--resultfld argument, it will default to the stem of the specified csv-filename. - The optional argument
--hrfallows you to specify a custom hrf saved as./code/hrf/CUSTOM-HRF.txt. This should be a single column text file with a sample frequency of 0.05s/sample. For more complicated FLOBS-style multi-function HRF's you can also specify multiple columns. If this argument is omitted, FSL's default double-gamma canonical HRF is used. This is often not a bad approximation, but a specific monkey HRF can be a little faster (there is code in Tracker-MRI to run an experiment that allows estimating the HRF on an individual basis). This argument default to FSL's standarddoublegamma, but we strongly recommend using a monkey HRF instead. - The optional argument
--MotionOutliersdetermines which motion outliers file is used. The default value issinglewhich means it uses the motion outlier file generated bybids_preprocessing_workflow.py. For a more conservative option you can specifymergedand it will pick the merged outliers files frombids_preprocessing_workflow.pywhich is a combination of outliers detected byrapidartandfsl_motion_correct. - The
--fwhmand--HighPassflags are optional to specify spatial smoothing (mm) and a high-pass filter (s). If not specified, they will be 2.0 mm and 50 s respectively.- example:
clear && ./code/bids_modelfit_level1_workflow.py --proj <ProjectName> --csv ./code/csv/ctcheckerboard.csv --contrasts ctcheckerboard |& tee ./logs/log-modelfit.txtorclear && ./code/bids_modelfit_level1_workflow.py --proj <ProjectName> --csv ./code/csv/checkerboard-ct-mapping.csv --contrasts ctcheckerboard --hrf ./code/hrf/custom-hrf.txt --resultfld my_unique_output --MotionOutliers single |& tee ./logs/log-modelfit.txt
- example:
- Log on to
lisa.surfsara.nl, go toNHP-BIDSdirectory, and runsbatch ./code/lisa/modelfit/<modelfit_script.sh>.
A level2 fixed effects analysis can be executed with ./code/bids_modelfit_level2_fixedfx_workflow.py. This workflow will take the results from level1 analysis and calculate the level2 statistics, which it will output to ./derivatives/modelfit/<contrast>/<registration_space>/level2/<resultfld>/fx/. Note that this analysis is specific for the contrast and registration space the level analysis was run with.
- The
--projflag is mandatory and indicates to which project these files belong. There should be a project folder inNHP-BIDS/projectswith the same label. This only hold for thenin-dataprotocolbranch for now but will become the default after testing. - The
--csvflag is mandatory and should point to the csv-file that is formatted as explained above. Here you can again specify a refsubject column in the csv. Make sure this is done correctly as it determines which brainmask file will be used in the analysis. - The
--contrastsflag is mandatory and depends on the experiment. Incode/contrasts/there are python modules for each set of contrasts. If you want to use the contrasts defined inctcheckerboard.py, for example, passctcheckerboardas the value for the--contrastsparameter. If you create your own contrasts, you just need your function to define the variable namedcontrasts. Since the code assumes that this will be a python module name, you cannot use dashes or spaces in the name. - The optional argument
--RegSpaceallows us to specify in which space the GLM is performed. The default isnmt, which means that the workflow will use the highpassed files that were aligned to theNMTtemplate space as input to the GLM. This puts all subjects in a common reference space, but it only works ifbids_warp2nmt_workflow.pyhas been executed already. The alternative optionnativewill use the preprocessed functional timeseries in their functional space. Be careful with this as there is no explicit check for mismatching functional spaces with multiple entries and statistics will be meaningless if not calculated in a common space for all functional runs. - The optional argument
--resultfldallows you to define the name of the working directory for this analysis. Re-running something with the same name uses existing intermediate results from working directories. Output directories are determined based on the contrast and registration space that are chosen allowing straightforward selection of results for level2 analysis. If you do not include a `--resultfld argument, it will default to the stem of the specified csv-filename. - The optional argument
--MotionOutliersdetermines which motion outliers level1 files are used. The default value isNonewhich means it ignores the whole mof selection and uses files in./derivatives/modelfit/<contrast>/<registration_space>/level1/selected/<sub>/<ses>/<run>. Note that such files are not automatically generated. You can manually select and copy or symlink as desired. If--MotionOutliersis specified assingleormergedit will take files from./derivatives/modelfit/<contrast>/<registration_space>/level1/mo-<motion_outlier_type>/<sub>/<ses/<run>instead.
In general, it is good practice to script your visualisation steps for reproducibility (and speed). We recommend doing this in project-specific github-repositories that house Jupyter notebooks for different visualisation tasks. If you ran the modelfit on warped data, the results will already be in the common template space and you can display it on both the NMT volume, and the corresponding surfaces. You can find all of these on the VS02 server. Alternatively you may want to display results on the individual anatomy and/or surfaces. In that case, you either use the reference-vols anatomicals, or have a look at our NHP-Freesurfer repository with Jupyter notebooks explain how to generate the surfaces from the individual monkey anatomical scans and project statistical result to the surface. For working in surface-space with Python, check out our adaptation of pyCortex called NHP-pyCortex.