Pipeline for processing functional data of the 7T MELAS study

public 1yr ago Version: 2 0 bookmarks

View Workflow

Help improve this workflow!

This workflow has been published but could be further improved with some additional meta data:

Keyword(s) in categories input, output, operation, topic

You can help improve this workflow by suggesting the addition or removal of keywords, suggest changes and report issues, or request to become a maintainer of the Workflow .

Snakemake workflow for processing functional data from the 7T MELAS study

Inputs:

participants.tsv with target subject IDs
For each target subject:
- Freesurfer processed data
- Functional data

Singularity containers required:

fMRIprep (for FreeSurfer, ADNI and FSL tools)
Connectome workbench

Authors

Roy AM Haast @royhaast

Usage

If you use this workflow in a paper, don't forget to give credits to the authors by citing the URL of this (original) repository and, if available, its DOI (see above).

Step 1: Obtain a copy of this workflow

Create a new github repository using this workflow as a template .
Clone the newly created repository to your local system, into the place where you want to perform the data analysis.

Step 2: Configure workflow

Configure the workflow according to your needs via editing the files in the config/ folder. Adjust config.yml to configure the workflow execution, and participants.tsv to specify your subjects.

Step 3: Install Snakemake

Install Snakemake using conda :

conda create -c bioconda -c conda-forge -n snakemake snakemake

For installation details, see the instructions in the Snakemake documentation .

Step 4: Execute workflow

Activate the conda environment:

conda activate snakemake

Test your configuration by performing a dry-run via

snakemake --use-singularity -n

Execute the workflow locally via

snakemake --use-singularity --cores $N

using $N cores or run it in a cluster environment via

snakemake --use-singularity --cluster qsub --jobs 100

snakemake --use-singularity --drmaa --jobs 100

If you are using Compute Canada, you can use the cc-slurm profile, which submits jobs and takes care of requesting the correct resources per job (including GPUs). Once it is set-up with cookiecutter, run:

snakemake --profile cc-slurm

Or, with neuroglia-helpers can get a 8-core, 32gb node and run locally there. First, get a node (default 8-core, 32gb, 3 hour limit):

regularInteractive

Then, run:

snakemake --use-singularity --cores 8 --resources mem=32000

See the Snakemake documentation for further details.

Step 5: Investigate results

After successful execution, you can create a self-contained interactive HTML report with all results via:

snakemake --report report.html

This report can, e.g., be forwarded to your collaborators. An example (using some trivial test data) can be seen here .

Step 6: Commit changes

Whenever you change something, don't forget to commit the changes back to your github copy of the repository:

git commit -a
git push

Step 7: Obtain updates from upstream

Whenever you want to synchronize your workflow copy with new developments from upstream, do the following.

Once, register the upstream repository in your local copy: git remote add -f upstream git@github.com:snakemake-workflows/{{cookiecutter.repo_name}}.git or git remote add -f upstream https://github.com/snakemake-workflows/{{cookiecutter.repo_name}}.git if you do not have setup ssh keys.
Update the upstream version: git fetch upstream .
Create a diff with the current version: git diff HEAD upstream/master workflow > upstream-changes.diff .
Investigate the changes: vim upstream-changes.diff .
Apply the modified diff via: git apply upstream-changes.diff .
Carefully check whether you need to update the config files: git diff HEAD upstream/master config . If so, do it manually, and only where necessary, since you would otherwise likely overwrite your settings and samples.

Step 8: Contribute back

In case you have also changed or added steps, please consider contributing them back to the original repository:

Fork the original repo to a personal or lab account.
Clone the fork to your local system, to a different place than where you ran your analysis.
Copy the modified files from your analysis to the clone of your fork, e.g., cp -r workflow path/to/fork . Make sure to not accidentally copy config file contents or sample sheets. Instead, manually update the example config files if necessary.
Commit and push your changes to your fork.
Create a pull request against the original repository.

Testing

TODO: create some test datasets

Code Snippets

import numpy as np
import pandas as pd
import nibabel as nib
from nibabel.nifti1 import Nifti1Image
from nilearn.input_data import NiftiLabelsMasker

epi = nib.load(snakemake.input.vol)
csf = nib.load(snakemake.input.csf)
wm = nib.load(snakemake.input.wm)

img = Nifti1Image(csf.get_fdata(),epi.affine,csf.header)
masker = NiftiLabelsMasker(labels_img=img, standardize=False)
time_series1 = masker.fit_transform(epi)

img = Nifti1Image(wm.get_fdata(),epi.affine,wm.header)
masker = NiftiLabelsMasker(labels_img=img, standardize=False)
time_series2 = masker.fit_transform(epi)

# Concatenate timeseries
df1 = pd.DataFrame({'CSF': time_series1[:,0],'WhiteMatter': time_series2[:,0]})

# Load movement parameters (and their derivatives)
names = ['X','Y','Z','RotX','RotY','RotZ']
df2 = pd.read_csv(snakemake.input.movreg, names=names, header=None, delim_whitespace=True)

# Put everything together (excluding derivatives) and write to .tsv file for 'ciftify_clean_img' 
df_concat = pd.concat([df2, df1], axis=1)
df_concat.to_csv(snakemake.output[0],sep='\t')

Python Pandas numpy nibabel nilearn From line 1 of scripts/extract_confounds.py

epi=$1
mni=$2

# Relevant transformations
gdc_warp=$3
mc_dir=$4
topup_dir=$5
bbr_mat=$6
fnirt_warp=$7

# Output paths
warped_mni=$8
out_dir=`dirname $warped_mni`

# Generate second warp: gradient distortion corrected to MNI
fslmaths $topup_dir/epi_fout.nii.gz -mul 0.0343097 $topup_dir/epi_fout_scaled.nii.gz
secondwarp=$out_dir/gdc_to_mni.nii.gz
convertwarp -s $topup_dir/epi_fout_scaled.nii.gz -d y- --premat=$bbr_mat --warp1=$fnirt_warp --ref=$mni --out=$secondwarp --relout --rel -v

# Generate seperate warp up to topup correction
intermediatewarp=$out_dir/gdc_to_topup.nii.gz
convertwarp -s $topup_dir/epi_fout_scaled.nii.gz -d y- --ref=${epi}_0000.nii.gz --out=$intermediatewarp --relout --rel -v

for file in $mc_dir/* ; do
    vol=`echo $file | rev | cut -d'_' -f1  | rev`
    epi_vol=${epi}_${vol}

    # Generate first warp: slicetiming corrected to gradient distortion corrected
    firstwarp=$out_dir/firstwarps/st_${vol}_to_gdc.nii.gz

    if [ ! -d $out_dir/firstwarps ] ; then mkdir -p $out_dir/firstwarps ; fi
    convertwarp --warp1=$gdc_warp --midmat=$file --ref=$epi_vol --out=$firstwarp --abs --relout -v

    # Combine with second warp: slicetiming corrected to MNI
    combinedwarp=$out_dir/combinedwarps/st_${vol}_to_MNI.nii.gz

    if [ ! -d $out_dir/combinedwarps ] ; then mkdir -p $out_dir/combinedwarps ; fi
    convertwarp --warp1=$firstwarp --warp2=$secondwarp --out=$combinedwarp --ref=$mni --rel --relout -v

    # Combine with intermediate warp: slicetiming corrected to topup correction
    topupwarp=$out_dir/topupwarps/st_${vol}_to_topup.nii.gz

    if [ ! -d $out_dir/topupwarps ] ; then mkdir -p $out_dir/topupwarps ; fi
    convertwarp --warp1=$firstwarp --warp2=$intermediatewarp --out=$topupwarp --ref=$epi_vol --rel --relout -v

    # Apply warps to MNI
    if [ ! -d $out_dir/warped_mni ] ; then mkdir -p $out_dir/warped_mni ; fi
    applywarp -i $epi_vol -r $mni -o $out_dir/warped_mni/epi_rest_${vol}_mni.nii.gz -w $combinedwarp --rel --interp=spline -v

    # Apply warps to topup correction
    if [ ! -d $out_dir/warped_topup ] ; then mkdir -p $out_dir/warped_topup ; fi    
    applywarp -i $epi_vol -r $epi_vol -o $out_dir/warped_topup/epi_rest_${vol}_topup.nii.gz -w $topupwarp --rel --interp=spline -v
done

# Merge warped output into 4D volume
fslmerge -tr $warped_mni `ls $out_dir/warped_mni/epi_rest_*_mni.nii.gz` 2.0
fslmerge -tr $out_dir/epi_rest_topup.nii.gz `ls $out_dir/warped_topup/epi_rest_*topup.nii.gz` 2.0

Shell From line 4 of scripts/onestepresampling.sh

shell:
    "{params.fs_setup} mri_convert {input.parc} {output.parc_anat} && "
    "applywarp -i {output.parc_anat} -r {params.mni} -w {input.warp} -o {output.parc_mni} -d int --interp=nn  && "
    "flirt -in {output.parc_anat} -ref {input.epi} -applyxfm -init {input.inverse_mat} -datatype int -interp nearestneighbour -out {output.parc_epi}"

SnakeMake From line 21 of rules/cifti.smk

shell:
    "{params.fs_setup} mri_binarize --i {input} --match 2 41 --erode 2 --o {output.wm} && "
    "mri_binarize --i {input} --match 4 43 --erode 1 --o {output.csf} && "
    "mri_binarize --i {input} --match 8 10 11 12 13 16 17 18 26 28 47 49 50 51 52 53 54 58 60 --o {output.rois} && "
    "mri_binarize --i {input} --match 0 --inv --o {output.brain} && "
    "fslmaths {output.rois} -mul {input} {output.rois}"

SnakeMake From line 40 of rules/cifti.smk

shell:
    "wb_command -volume-label-import {input.rois} {input.labels} {output}"

SnakeMake From line 54 of rules/cifti.smk

shell:
    "wb_command -cifti-create-label {output} -volume {input.rois} {input.rois} -left-label {input.lh_mmp} -right-label {input.rh_mmp}"

SnakeMake From line 65 of rules/cifti.smk

shell:
    "wb_command -cifti-create-dense-timeseries {output} -volume {input.vol} {input.rois} -left-metric {input.lh} -right-metric {input.rh} -timestep 2.0"

SnakeMake From line 80 of rules/cifti.smk

shell:
    "wb_command -cifti-create-dense-timeseries {output} -volume {input.vol} {input.rois} -left-metric {input.lh} -right-metric {input.rh} -timestep 2.0"

SnakeMake From line 95 of rules/cifti.smk

script:
    "../../scripts/extract_confounds.py"

SnakeMake From line 110 of rules/cifti.smk

shell:
    "ciftify_clean_img --verbose --output-file={output} --confounds-tsv={input.confounds} "
    "--smooth-fwhm=3 --left-surface {input.lh_surf} --right-surface {input.rh_surf} {params} {input.dtseries} &> {log}"

SnakeMake From line 129 of rules/cifti.smk

shell:
    'ciftify_clean_img --verbose --output-file={output} --confounds-tsv={input.confounds} {params} {input.dtseries} &> {log}'

SnakeMake From line 145 of rules/cifti.smk

shell:
    "wb_command -cifti-parcellate {input.dtseries} {input.dlabels} COLUMN {output}"

SnakeMake From line 158 of rules/cifti.smk

shell:
    "wb_command -cifti-correlation {input} {output}"

SnakeMake From line 169 of rules/cifti.smk

shell:
    "{params.fs_setup} bbregister --s {wildcards.subject} --mov {input.epi} --reg {output.reg} --o {output.vol} {params.optional} && "
    "tkregister2 --noedit --reg {output.reg} --mov {input.epi} --targ {input.t1w} --fslregout {output.mat} &> {log}"      

SnakeMake From line 21 of rules/coregistration.smk

shell:
    "convert_xfm -omat {output.inverse_mat} -inverse {input.mat} && "
    "flirt -interp spline -in {input.t1w} -ref {input.epi} -applyxfm -init {output.inverse_mat} -out {output.t1w_coreg}"

SnakeMake From line 38 of rules/coregistration.smk

shell:
    "{params.fs_setup} mri_vol2vol --mov {input.epi} --targ {input.targ} --o {output} --reg {input.bbr} --no-resample"

SnakeMake From line 55 of rules/coregistration.smk

shell:
    "flirt {params.optional} -in {input} -ref {params.MNI} -omat {output} -out {params.prefix} &> {log}"

SnakeMake From line 74 of rules/coregistration.smk

shell:
    "fnirt --in={input.t1w} --ref={params.MNI} --aff={input.affine} --fout={output.fout} --jout={output.jout} "
    "--refout={output.refout} --iout={output.iout} --logout={output.logout} --intout={output.intout} --cout={output.cout} -v &> {log}"

SnakeMake From line 98 of rules/coregistration.smk

shell:
    "invwarp -w {input.warp} -o {output} -r {input.t1w} --noconstraint"

SnakeMake From line 113 of rules/coregistration.smk

shell:
    "bash scripts/onestepresampling.sh {params.epi} {params.mni} {input.gdc} {params.mat} "
    "{input.topup} {input.bbr} {input.fnirt} {output.warped_mni} &> {log}"

SnakeMake From line 136 of rules/coregistration.smk

shell:
    "python /opt/ICA-AROMA/ICA_AROMA.py -o `dirname {output.den}` -i {input.vol} -mc {input.mc} "
    "-a {input.bbr} -w {input.fnirt} -m {input.mask} -tr 2 -den nonaggr -overwrite"

SnakeMake From line 154 of rules/coregistration.smk

shell:
    "applywarp -i {input.vol} -r {input.mni} -o {output} -w {input.warp} --rel --interp=spline -v"

SnakeMake From line 168 of rules/coregistration.smk

shell:
    "TR=`fslval {input[0]} pixdim4` && "
    "3dTshift -prefix {output.corrected} -tpattern @{params.slicetiming} -TR $TR -tzero 0.0 {input} && "
    "fslsplit {output.corrected} `echo {output.splitted} | rev | cut -d'_' -f2-  | rev`_ -t &> {log}"

SnakeMake From line 15 of rules/preprocessing.smk

shell:
    "fslroi {input} {output.firstvol} 0 1 && "
    "procGradCorrect -i {output.firstvol} -g {params.coeff} -w {output.warp} -s {params.scratch} &> {log}"

SnakeMake From line 35 of rules/preprocessing.smk

shell:
    "applywarp -i {input.epi} -o {output.unwarped} -r {input.epi} -w {input.warp} --abs --interp=spline && "
    "fslroi {output.unwarped} {output.jacobian} 0 1 && "
    "reg_jacobian -ref {output.jacobian} -def {input.warp} -jac {output.jacobian} && "
    "fslmaths {output.jacobian} -mul -1 -abs {output.jacobian} && "
    "fslmaths {output.unwarped} -mul {output.jacobian} {output.corrected} && "
    "fslcpgeom {input.epi} {output.corrected}"

SnakeMake From line 53 of rules/preprocessing.smk

shell:
    "fslroi {input[0]} {output.first_vols} 0 5 && "
    "fslmerge -tr {output.concat_out} {output.first_vols} {input[1]} 2"

SnakeMake From line 72 of rules/preprocessing.smk

shell:
    "mkdir -p {output} && "
    "topup --imain={input} --datain={params.acquisition} --config={params.topup} "
    "--out={output}/epi --fout={output}/epi_fout --iout={output}/epi_iout -v &> {log}"

SnakeMake From line 89 of rules/preprocessing.smk

shell:
    "mcflirt -in {input} -out {output.vol} {params.optional} &> {log}"

SnakeMake From line 108 of rules/preprocessing.smk

shell:
    "fslroi {input.epi} {output.firstvol} 0 1 && "
    "applytopup --imain={output.firstvol} --inindex=1 --datain={params.acquisition} --topup={input.topup}/epi --out={output.firstvol_gdc} --method=jac -v &> {log} && "
    "applytopup --imain={input.epi} --inindex=1 --datain={params.acquisition} --topup={input.topup}/epi --out={output.gdc} --method=jac -v &> {log}"

SnakeMake From line 128 of rules/preprocessing.smk

shell:
    "bet {input.mask} {output} -f 0.3 -n -m -R"

SnakeMake From line 137 of rules/preprocessing.smk

shell:
    "{params.fs_setup} mris_convert {input} {output}"

SnakeMake From line 10 of rules/surface.smk

shell:
    "wb_command -surface-cortex-layer {input.white} {input.pial} 0.5 {output}"      

SnakeMake From line 20 of rules/surface.smk

shell: 
    "{params.fs_setup} mri_info {input} --tkr2scanner > {output}"

SnakeMake From line 32 of rules/surface.smk

shell:
    "wb_command -convert-affine -from-flirt {input.bbr} {input.t1w} {input.epi} -to-world {output}"

SnakeMake From line 44 of rules/surface.smk

shell: 
    "wb_command -surface-apply-affine {input.surf} {input.tkr2scanner} {output}"    

SnakeMake From line 55 of rules/surface.smk

shell:     
    "wb_command -surface-apply-affine {input.surf} {input.scanner2epi} {output}"

SnakeMake From line 66 of rules/surface.smk

shell:
    "wb_command -surface-apply-warpfield {input.surf} {input.inverse_warp} {output} -fnirt {input.warp}"

SnakeMake From line 78 of rules/surface.smk

shell:    
    "wb_command -surface-resample {input.surf} {input.sphere_old} {input.sphere_new} BARYCENTRIC {output}"

SnakeMake From line 91 of rules/surface.smk

shell:     
    "wb_command -surface-generate-inflated {input} {output.inflated} {output.very_inflated} {params}"   

SnakeMake From line 103 of rules/surface.smk

shell: 
    "wb_command -volume-to-surface-mapping {input.epi} {input.surf} {output} -trilinear"

SnakeMake From line 119 of rules/surface.smk

shell: 
    "wb_command -volume-to-surface-mapping {input.epi} {input.surf} {output} -trilinear"

SnakeMake From line 133 of rules/surface.smk

shell:     
    "wb_command -metric-resample {input.metric} {input.sphere_old} {input.sphere_new} ADAP_BARY_AREA {output} -area-surfs {input.area_old} {input.area_new}"

SnakeMake From line 149 of rules/surface.smk

shell:     
    "wb_command -metric-resample {input.metric} {input.sphere_old} {input.sphere_new} ADAP_BARY_AREA {output} -area-surfs {input.area_old} {input.area_new}"