intial commit

.gitattributes

@@ -2,4 +2,6 @@
 * annex.backend=MD5E
 **/.git* annex.largefiles=nothing
 CHANGELOG.md annex.largefiles=nothing
-README.md annex.largefiles=nothing
+README.md annex.largefiles=nothing
+requirements.txt annex.largefiles=nothing
+highspeed-masks.Rproj annex.largefiles=nothing

.gitignore

@@ -0,0 +1,5 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+code/*/*.html

code/docs/highspeed-masks-docs.Rmd

@@ -0,0 +1,49 @@
+---
+title: "Faster than thought: Detecting sub-second activation sequences with sequential fMRI pattern analysis"
+subtitle: "Short project title: highspeed"
+author:
+  - Lennart Wittkuhn^[Max Planck Institute for Human Development, wittkuhn@mpib-berlin.mpg.de]
+  - Nicolas W. Schuck^[Max Planck Institute for Human Development, schuck@mpib-berlin.mpg.de]
+date: "Last update: `r format(Sys.time(), '%d %B, %Y')`"
+output:
+ html_document:
+    toc: true
+    self_contained: true
+    toc_float: true
+    toc_depth: 3
+    number_sections: true
+    highlight: pygments
+    theme: cosmo
+    df_print: paged
+    fig_caption: true
+fig.align: "center"
+header-includes:
+  - \usepackage{fontspec}
+  - \setmainfont{AgfaRotisSansSerif}
+email: wittkuhn@mpib-berlin.mpg.de
+---
+
+```{r, echo=FALSE, message=FALSE, include=FALSE}
+if (!requireNamespace("pacman")) install.packages("pacman")
+packages_cran <- c("here")
+pacman::p_load(char = packages_cran)
+if (basename(here::here()) == "highspeed"){
+  path_root = here::here("highspeed-masks")
+} else {
+  path_root = here::here()
+}
+```
+
+# Creating binary masks
+
+### Code and software
+
+#### Creating binary masks using `highspeed-masks.py`
+
+```{python, echo=TRUE, code=readLines(file.path(path_root, "code", "masks", "highspeed-masks.py")), eval=FALSE, python.reticulate=FALSE}
+```
+
+#### Plotting masked data using `highspeed-masks-plot.py`
+
+```{python, echo=TRUE, code=readLines(file.path(path_root, "code", "masks", "highspeed-masks-plot.py")), eval=FALSE, python.reticulate=FALSE}
+```

code/masks/highspeed-masks-plot.py

@@ -0,0 +1,91 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# ======================================================================
+# SCRIPT INFORMATION:
+# ======================================================================
+# SCRIPT: PLOT MASKS FOR ILLUSTRATION
+# PROJECT: HIGHSPEED
+# WRITTEN BY LENNART WITTKUHN, 2020
+# CONTACT: WITTKUHN AT MPIB HYPHEN BERLIN DOT MPG DOT DE
+# MAX PLANCK RESEARCH GROUP NEUROCODE
+# MAX PLANCK INSTITUTE FOR HUMAN DEVELOPMENT
+# MAX PLANCK UCL CENTRE FOR COMPUTATIONAL PSYCHIATRY AND AGEING RESEARCH
+# LENTZEALLEE 94, 14195 BERLIN, GERMANY
+# ======================================================================
+# IMPORT RELEVANT PACKAGES
+# ======================================================================
+from nilearn import plotting
+from nilearn import image
+import os
+from os.path import join as opj
+import glob
+import re
+import matplotlib.pyplot as plt
+import statistics
+
+sub = 'sub-01'
+path_root = os.getcwd()
+path_anat = opj(path_root, 'data', 'bids', sub, '*', 'anat', '*.nii.gz')
+anat = sorted(glob.glob(path_anat))[0]
+path_patterns = opj(path_root, 'data', 'decoding', sub, 'data', '*.nii.gz')
+path_patterns = sorted(glob.glob(path_patterns))
+path_patterns = [x for x in path_patterns if 'hpc' not in x]
+path_union = opj(path_root, 'data', 'decoding', sub, 'masks', '*union.nii.gz')
+path_union = glob.glob(path_union)
+
+# path_fmriprep = opj(path_root, 'derivatives', 'fmriprep')
+# path_masks = opj(path_root, 'derivatives', 'decoding', sub)
+# path_anat = opj(path_fmriprep, sub, 'anat', sub + '_desc-preproc_T1w.nii.gz')
+# path_visual = opj(path_masks, 'masks', '*', '*.nii.gz')
+# vis_mask = glob.glob(path_visual)[0]
+# vis_mask_smooth = image.smooth_img(vis_mask, 4)
+
+plotting.plot_roi(path_union[0], bg_img=anat,
+                  cut_coords=[30, 10, 15],
+                  title="Region-of-interest", black_bg=True,
+                  display_mode='ortho', cmap='red_transparent',
+                  draw_cross=False)
+
+# calculate average patterns across all trials:
+mean_patterns = [image.mean_img(i) for i in path_patterns]
+# check the shape of the mean patterns (should be 3D):
+[print(image.get_data(i).shape) for i in mean_patterns]
+# extract labels of patterns
+labels = [re.search('union_(.+?).nii.gz', i).group(1) for i in path_patterns]
+
+
+# function used to plot individual patterns:
+def plot_patterns(pattern, name):
+    display = plotting.plot_stat_map(
+            pattern,
+            bg_img=anat,
+            #cut_coords=[30, 29, -6],
+            title=name,
+            black_bg=True,
+            colorbar=True,
+            display_mode='ortho',
+            draw_cross=False
+            )
+    path_save = opj(path_root, 'figures', 'pattern_{}.pdf').format(name)
+    display.savefig(filename=path_save)
+    display.close()
+    return(display)
+
+
+# plot individual patterns and save coordinates:
+coords = []
+for pattern, name in zip(mean_patterns, labels):
+    display = plot_patterns(pattern, name)
+    coords.append(display.cut_coords)
+# mean_coords = [sum(x)/len(x) for x in zip(*coords)]
+mean_coords = [statistics.median(x) for x in zip(*coords)]
+
+# create subplot with all patterns using mean coordinates:
+fig, axes = plt.subplots(nrows=len(path_patterns), ncols=1, figsize=(14, 20))
+for pattern, name, ax in zip(mean_patterns, labels,  axes):
+    display = plotting.plot_stat_map(
+            pattern, bg_img=anat, cut_coords=mean_coords, title=name,
+            black_bg=True, colorbar=True, display_mode='ortho',
+            draw_cross=False, axes=ax, symmetric_cbar=True, vmax=1)
+fig.subplots_adjust(wspace=0, hspace=0)
+fig.savefig(opj(path_root, 'figures', 'pattern_all.pdf'), bbox_inches='tight')

code/masks/highspeed-masks.py

@@ -0,0 +1,235 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# ======================================================================
+# SCRIPT INFORMATION:
+# ======================================================================
+# SCRIPT: CREATE BINARIZED MASKS FROM SEGMENTED FUNCTIONAL IMAGES
+# PROJECT: HIGHSPEED
+# WRITTEN BY LENNART WITTKUHN, 2018 - 2020
+# CONTACT: WITTKUHN AT MPIB HYPHEN BERLIN DOT MPG DOT DE
+# MAX PLANCK RESEARCH GROUP NEUROCODE
+# MAX PLANCK INSTITUTE FOR HUMAN DEVELOPMENT
+# MAX PLANCK UCL CENTRE FOR COMPUTATIONAL PSYCHIATRY AND AGEING RESEARCH
+# LENTZEALLEE 94, 14195 BERLIN, GERMANY
+# ======================================================================
+# IMPORT RELEVANT PACKAGES
+# ======================================================================
+# import basic libraries:
+import os
+import sys
+import warnings
+from os.path import join as opj
+# import nipype libraries:
+from nipype.interfaces.utility import IdentityInterface
+from nipype.interfaces.io import SelectFiles, DataSink
+from nipype.pipeline.engine import Workflow, Node, MapNode
+# import libraries for bids interaction:
+from bids.layout import BIDSLayout
+# import freesurfer interfaces:
+from nipype.interfaces.freesurfer import Binarize
+# import fsl interfaces:
+from niflow.nipype1.workflows.fmri.fsl import create_susan_smooth
+# ======================================================================
+# ENVIRONMENT SETTINGS (DEALING WITH ERRORS AND WARNINGS):
+# ======================================================================
+# set the fsl output type environment variable:
+os.environ['FSLOUTPUTTYPE'] = 'NIFTI_GZ'
+# set the freesurfer subject directory:
+os.environ['SUBJECTS_DIR'] = '/home/mpib/wittkuhn/tools/freesurfer/'
+# deal with nipype-related warnings:
+os.environ['TF_CPP_MIN_LOG_LEVEL'] = "3"
+# filter out warnings related to the numpy package:
+warnings.filterwarnings("ignore", message="numpy.dtype size changed*")
+warnings.filterwarnings("ignore", message="numpy.ufunc size changed*")
+# ======================================================================
+# DEFINE PBS CLUSTER JOB TEMPLATE (NEEDED WHEN RUNNING ON THE CLUSTER):
+# ======================================================================
+job_template = """
+#PBS -l walltime=5:00:00
+#PBS -j oe
+#PBS -o /home/mpib/wittkuhn/highspeed/logs/masks
+#PBS -m n
+#PBS -v FSLOUTPUTTYPE=NIFTI_GZ
+source /etc/bash_completion.d/virtualenvwrapper
+workon glm
+module load fsl/5.0
+module load freesurfer/6.0.0
+"""
+# ======================================================================
+# DEFINE PATHS AND SUBJECTS
+# ======================================================================
+# define paths depending on the operating system (OS) platform:
+path_root = None
+sub_list = None
+if 'darwin' in sys.platform:
+    path_root = opj('/Users', 'wittkuhn', 'Volumes', 'tardis', 'highspeed')
+    sub_list = ['sub-01']
+elif 'linux' in sys.platform:
+    path_root = opj('/home', 'mpib', 'wittkuhn', 'highspeed')
+    # grab the list of subjects from the bids data set:
+    layout = BIDSLayout(opj(path_root, 'bids'))
+    # get all subject ids:
+    sub_list = sorted(layout.get_subjects())
+    # create a template to add the "sub-" prefix to the ids
+    sub_template = ['sub-'] * len(sub_list)
+    # add the prefix to all ids:
+    sub_list = ["%s%s" % t for t in zip(sub_template, sub_list)]
+    # if user defined to run specific subject
+    sub_list = sub_list[int(sys.argv[1]):int(sys.argv[2])]
+# ======================================================================
+# DEFINE NODE: INFOSOURCE
+# ======================================================================
+# define the infosource node that collects the data:
+infosource = Node(IdentityInterface(
+    fields=['subject_id']), name='infosource')
+# let the node iterate (parallelize) over all subjects:
+infosource.iterables = [('subject_id', sub_list)]
+# ======================================================================
+# DEFINE SELECTFILES NODE
+# ======================================================================
+# define all relevant files paths:
+templates = dict(
+    func=opj(path_root, 'derivatives', 'fmriprep', '{subject_id}', '*',
+             'func', '*space-T1w*preproc_bold.nii.gz'),
+    func_parc=opj(path_root, 'derivatives', 'fmriprep', '{subject_id}',
+                  '*', 'func', '*space-T1w*aparcaseg_dseg.nii.gz'),
+    wholemask=opj(path_root, 'derivatives', 'fmriprep', '{subject_id}',
+                  '*', 'func', '*space-T1w*brain_mask.nii.gz'),
+)
+# define the selectfiles node:
+selectfiles = Node(SelectFiles(templates, sort_filelist=True),
+                   name='selectfiles')
+# set expected thread and memory usage for the node:
+selectfiles.interface.num_threads = 1
+selectfiles.interface.estimated_memory_gb = 0.1
+# selectfiles.inputs.subject_id = 'sub-20'
+# selectfiles_results = selectfiles.run()
+# ======================================================================
+# DEFINE CREATE_SUSAN_SMOOTH WORKFLOW NODE
+# ======================================================================
+# define the susan smoothing node and specify the smoothing fwhm:
+susan = create_susan_smooth()
+# set the smoothing kernel:
+susan.inputs.inputnode.fwhm = 4
+# set expected thread and memory usage for the nodes:
+susan.get_node('inputnode').interface.num_threads = 1
+susan.get_node('inputnode').interface.estimated_memory_gb = 0.1
+susan.get_node('median').interface.num_threads = 1
+susan.get_node('median').interface.estimated_memory_gb = 3
+susan.get_node('mask').interface.num_threads = 1
+susan.get_node('mask').interface.estimated_memory_gb = 3
+susan.get_node('meanfunc2').interface.num_threads = 1
+susan.get_node('meanfunc2').interface.estimated_memory_gb = 3
+susan.get_node('merge').interface.num_threads = 1
+susan.get_node('merge').interface.estimated_memory_gb = 3
+susan.get_node('multi_inputs').interface.num_threads = 1
+susan.get_node('multi_inputs').interface.estimated_memory_gb = 3
+susan.get_node('smooth').interface.num_threads = 1
+susan.get_node('smooth').interface.estimated_memory_gb = 3
+susan.get_node('outputnode').interface.num_threads = 1
+susan.get_node('outputnode').interface.estimated_memory_gb = 0.1
+# ======================================================================
+# DEFINE BINARIZE NODE
+# ======================================================================
+mask_visual_labels = [
+    1005, 2005,  # cuneus
+    1011, 2011,  # lateral occipital
+    1021, 2021,  # pericalcarine
+    1029, 2029,  # superioparietal
+    1013, 2013,  # lingual
+    1008, 2008,  # inferioparietal
+    1007, 2007,  # fusiform
+    1009, 2009,  # inferiotemporal
+    1016, 2016,  # parahippocampal
+    1015, 2015,  # middle temporal
+]
+mask_hippocampus_labels = [
+    17, 53,  # left and right hippocampus
+]
+mask_mtl_labels = [
+    17, 53,  # left and right hippocampus
+    1016, 2016,  # parahippocampal
+    1006, 2006,  # ctx-entorhinal
+]
+# function: use freesurfer mri_binarize to threshold an input volume
+mask_visual = MapNode(
+        interface=Binarize(), name='mask_visual', iterfield=['in_file'])
+# input: match instead of threshold
+mask_visual.inputs.match = mask_visual_labels
+# optimize the efficiency of the node:
+mask_visual.plugin_args = {'qsub_args': '-l nodes=1:ppn=1', 'overwrite': True}
+mask_visual.plugin_args = {'qsub_args': '-l mem=100MB', 'overwrite': True}
+# function: use freesurfer mri_binarize to threshold an input volume
+mask_hippocampus = MapNode(
+        interface=Binarize(), name='mask_hippocampus', iterfield=['in_file'])
+# input: match instead of threshold
+mask_hippocampus.inputs.match = mask_hippocampus_labels
+# optimize the efficiency of the node:
+mask_hippocampus.plugin_args = {
+    'qsub_args': '-l nodes=1:ppn=1', 'overwrite': True}
+mask_hippocampus.plugin_args = {
+    'qsub_args': '-l mem=100MB', 'overwrite': True}
+# function: use freesurfer mri_binarize to threshold an input volume
+mask_mtl = MapNode(
+        interface=Binarize(), name='mask_mtl', iterfield=['in_file'])
+# input: match instead of threshold
+mask_mtl.inputs.match = mask_mtl_labels
+# optimize the efficiency of the node:
+mask_mtl.plugin_args = {'qsub_args': '-l nodes=1:ppn=1', 'overwrite': True}
+mask_mtl.plugin_args = {'qsub_args': '-l mem=100MB', 'overwrite': True}
+# ======================================================================
+# CREATE DATASINK NODE
+# ======================================================================
+# create a node of the function:
+datasink = Node(DataSink(), name='datasink')
+# assign the path to the base directory:
+datasink.inputs.base_directory = opj(path_root, 'derivatives', 'masks')
+# create a list of substitutions to adjust the filepaths of datasink:
+substitutions = [('_subject_id_', '')]
+# assign the substitutions to the datasink command:
+datasink.inputs.substitutions = substitutions
+# determine whether to store output in parameterized form:
+datasink.inputs.parameterization = True
+# set expected thread and memory usage for the node:
+datasink.interface.num_threads = 1
+datasink.interface.estimated_memory_gb = 0.2
+# ======================================================================
+# DEFINE WORKFLOW PIPELINE:
+# ======================================================================
+# initiation of the 1st-level analysis workflow:
+wf = Workflow(name='masks')
+# stop execution of the workflow if an error is encountered:
+wf.config = {'execution': {'stop_on_first_crash': True}}
+# define the base directory of the workflow:
+wf.base_dir = opj(path_root, 'work')
+# connect infosource to selectfiles node:
+wf.connect(infosource, 'subject_id', selectfiles, 'subject_id')
+# connect functional files to smoothing workflow:
+wf.connect(selectfiles, 'func', susan, 'inputnode.in_files')
+wf.connect(selectfiles, 'wholemask', susan, 'inputnode.mask_file')
+wf.connect(susan, 'outputnode.smoothed_files', datasink, 'smooth')
+# connect segmented functional files to visual mask node
+wf.connect(selectfiles, 'func_parc', mask_visual, 'in_file')
+wf.connect(mask_visual, 'binary_file', datasink, 'mask_visual.@binary')
+wf.connect(mask_visual, 'count_file', datasink, 'mask_visual.@count')
+# connect segmented functional files to hippocampus node
+wf.connect(selectfiles, 'func_parc', mask_hippocampus, 'in_file')
+wf.connect(mask_hippocampus, 'binary_file', datasink, 'mask_hippocampus.@binary')
+wf.connect(mask_hippocampus, 'count_file', datasink, 'mask_hippocampus.@count')
+# connect segmented functional files to mtl node
+wf.connect(selectfiles, 'func_parc', mask_mtl, 'in_file')
+wf.connect(mask_mtl, 'binary_file', datasink, 'mask_mtl.@binary')
+wf.connect(mask_mtl, 'count_file', datasink, 'mask_mtl.@count')
+# ======================================================================
+# WRITE GRAPH AND EXECUTE THE WORKFLOW
+# ======================================================================
+# write the graph:
+wf.write_graph(graph2use='colored', simple_form=True)
+# set the maximum resources the workflow can utilize:
+# args_dict = {'status_callback' : log_nodes_cb}
+# execute the workflow depending on the operating system:
+if 'darwin' in sys.platform:
+    # will execute the workflow using all available cpus:
+    wf.run(plugin='MultiProc')
+elif 'linux' in sys.platform:
+    wf.run(plugin='PBS', plugin_args=dict(template=job_template))

highspeed-masks.Rproj

@@ -0,0 +1,13 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX

requirements.txt

@@ -0,0 +1,42 @@
+bids-validator==1.3.12
+certifi==2019.11.28
+chardet==3.0.4
+Click==7.0
+decorator==4.4.1
+docopt==0.6.2
+etelemetry==0.1.2
+filelock==3.0.12
+future==0.18.2
+idna==2.8
+isodate==0.6.0
+joblib==0.14.1
+lxml==4.4.2
+networkx==2.4
+neurdflib==5.0.1
+nibabel==3.0.0
+niflow-nipype1-workflows==0.0.4
+nilearn==0.6.0
+nipype==1.4.0
+num2words==0.5.10
+numpy==1.18.0
+packaging==20.0
+pandas==0.25.3
+patsy==0.5.1
+pkg-resources==0.0.0
+prov==1.5.3
+pybids==0.10.0
+pydot==1.4.1
+pydotplus==2.0.2
+pyparsing==2.4.6
+python-dateutil==2.8.1
+pytz==2019.3
+rdflib==4.2.2
+requests==2.22.0
+scikit-learn==0.22.1
+scipy==1.4.1
+simplejson==3.17.0
+six==1.13.0
+sklearn==0.0
+SQLAlchemy==1.3.12
+traits==5.2.0
+urllib3==1.25.7