diff --git a/dbsc/config/snakebids.yml b/dbsc/config/snakebids.yml index 240093e5..c339875e 100644 --- a/dbsc/config/snakebids.yml +++ b/dbsc/config/snakebids.yml @@ -3,48 +3,54 @@ output_dir: "/path/to/output/dir" # don't use "." snakebids_dir: "." # DO NOT MODIFY - used to refer to files need by workflow running with CLI -debug: False # Enable printing of debug statements during parsing -- disable if generating dag visualizations +debug: False # Enable printing of debug statements during parsing -- disable if generating dag visualizations derivatives: False # Search in bids/derivatives if True; can also be path(s) to derivatives datasets # List of analysis levels in bids app analysis_levels: &analysis_levels - participant - - group # Mapping from analysis_level to set of target rules or files targets_by_analysis_level: participant: - "all" # if "", then the first rule is run - group: - - "all_group_tsv" # this configures the pybids grabber - create an entry for each type of input you want to grab # indexed by name of input dictionary for each input is passed directly to pybids get() # https://bids-standard.github.io/pybids/generated/bids.layout.BIDSLayout.html#bids.layout.BIDSLayout.get pybids_inputs: - T1w: + dwi: filters: - suffix: "T1w" + suffix: "dwi" extension: ".nii.gz" - datatype: "anat" - space: null + datatype: "dwi" + space: "T1w" + wildcards: + - subject + - session + mask: + filters: + suffix: "mask" + extension: ".nii.gz" + datatype: "dwi" + space: "T1w" wildcards: - subject - session # Configuration for the command-line parameters to make available # passed on the argparse add_argument() -parse_args: +parse_args: -#--- core BIDS-app options --- (do not modify below) ---# +#--- core BIDS-app options --- (do not modify below) ---# bids_dir: - help: The directory with the input dataset formatted according to the + help: The directory with the input dataset formatted according to the BIDS standard. output_dir: - help: The directory where the output files should be stored. If you are - running group level analysis, this folder should be prepopulated + help: The directory where the output files should be stored. If you are + running group level analysis, this folder should be prepopulated with the results of the participant level analysis. analysis_level: @@ -62,59 +68,64 @@ parse_args: --exclude_participant_label: help: The label(s) of the participant(s) that should be excluded. The label corresponds to sub- from the BIDS spec (so it does - not include "sub-"). If this parameter is not provided, all subjects + not include "sub-"). If this parameter is not provided, all subjects will be analyzed. Multiple participants can be specified with a space sepearated list. nargs: "+" -#-----------------------------------------------------# +#-----------------------------------------------------# #--- additional BIDS-app options --- (add in below) --# --shells: - help: '(Mrtrix3) specify one or more b-values to use during processing, as - a comma-separated list of the desired approximate b-values (b-values - are clustered to allow for small deviations). Note that some - commands are incompatible with multiple b-values, and will report an - error if more than one b-value is provided. + help: '(Mrtrix3) specify one or more b-values to use during processing, as + a comma-separated list of the desired approximate b-values (b-values + are clustered to allow for small deviations). Note that some + commands are incompatible with multiple b-values, and will report an + error if more than one b-value is provided. WARNING: note that, even though b=0 volumes are never referred to as - shells in the literature, they still have to be explicitly included - in the list of b-values as provided to the -shell option! Several - algorithms which include the b=0 volumes in their computations may + shells in the literature, they still have to be explicitly included + in the list of b-values as provided to the -shell option! Several + algorithms which include the b=0 volumes in their computations may otherwise return an undesired reuslt.' nargs: "+" --lmax: - help: '(Mrtrix3) the maximum spherical harmonic order for the output - FOD(s). For algorithms with multiple outputs, this should be - provided as a comma-sperated list of integers, one for each output - image; for single-output algorithms, only a single integer - should be provided. If omitted, commands will use the lmax of the - corresponding response function (i.e. based on its number of + help: '(Mrtrix3) the maximum spherical harmonic order for the output + FOD(s). For algorithms with multiple outputs, this should be + provided as a comma-sperated list of integers, one for each output + image; for single-output algorithms, only a single integer + should be provided. If omitted, commands will use the lmax of the + corresponding response function (i.e. based on its number of coefficients), up to a maximum of 8.' nargs: "+" --step: help: '(Mrtrix3) set the step size of the algorithm in mm. Should set to - 4 steps per voxel (e.g. 1/4 * voxel size) in order to sample more + 4 steps per voxel (e.g. 1/4 * voxel size) in order to sample more frequently in compact region. (default: 0.35mm)' default: 0.35 --sl_count: - help: '(Mrtrix3) set the desired number of streamlines to be selected, - after all selection criteria has been applied (i.e. - inclusion/exclusion ROIS, min/max length, etc. Streamlines will be - seeded until this number of streamlines have been selected, or the - maximum allowed number of seeds has been exceeded. Set to zero to + help: '(Mrtrix3) set the desired number of streamlines to be selected, + after all selection criteria has been applied (i.e. + inclusion/exclusion ROIS, min/max length, etc. Streamlines will be + seeded until this number of streamlines have been selected, or the + maximum allowed number of seeds has been exceeded. Set to zero to disable, which will result in streamlines being seeded until - the maximum allowed number of seeds has been reached. + the maximum allowed number of seeds has been reached. (default: 20,000,000 streamlines)' default: 20000000 --radial_search: - help: '(Mrtrix3) perform a radial search from each streamline endpoint to - locate the nearest node. Argument is the maximum radius in mm; if no - node is found within this radius, the streamline endpoint is not + help: '(Mrtrix3) perform a radial search from each streamline endpoint to + locate the nearest node. Argument is the maximum radius in mm; if no + node is found within this radius, the streamline endpoint is not assigned to any node. (default: 1.5mm)' default: 1.5 + + --responsemean_dir: + help: 'Provide directory containing average response functions. If not + provided, one will be computed from the subjects in the input + directory.' #-----------------------------------------------------# # Workflow specific config @@ -140,6 +151,6 @@ freesurfer: singularity: freesurfer: "docker://pwighton/freesurfer:7.2.0" neuroglia-core: "docker://khanlab/neuroglia-core:latest" - mrtpipelines: "docker://kaitj/mrtpipelines:0.1.6" + mrtrix: "docker://brainlife/mrtrix:3.0.3" -fs_license: /path/to/fs/license \ No newline at end of file +fs_license: /path/to/fs/license diff --git a/dbsc/run.py b/dbsc/run.py index ebaced76..7d23de77 100644 --- a/dbsc/run.py +++ b/dbsc/run.py @@ -1,9 +1,11 @@ import os + from snakebids.app import SnakeBidsApp + def get_parser(): """Exposes parser for sphinx doc generation, cwd is the docs dir""" - app = SnakeBidsApp('../surftexture',skip_parse_args=True) + app = SnakeBidsApp("../dbsc", skip_parse_args=True) return app.parser @@ -13,4 +15,4 @@ def main(): if __name__ == "__main__": - main() \ No newline at end of file + main() diff --git a/dbsc/workflow/Snakefile b/dbsc/workflow/Snakefile index 542d624c..db973afc 100644 --- a/dbsc/workflow/Snakefile +++ b/dbsc/workflow/Snakefile @@ -12,15 +12,18 @@ config.update( pybids_inputs=config["pybids_inputs"], derivatives=config["derivatives"], participant_label=config["participant_label"], + exclude_participant_label=config["exclude_participant_label"], + use_bids_inputs=False, ) ) # this adds constraints to the bids naming -wildcard_constraints: **snakebids.get_wildcard_constraints(\ - config["pybids_inputs"]\ -) +wildcard_constraints: + **snakebids.get_wildcard_constraints(config["pybids_inputs"]) + #---- end snakebids boilerplate ------------------------------------------------ # Rules include: "rules/freesurfer.smk" -include: "rules/zona_bb_subcortex" \ No newline at end of file +include: "rules/zona_bb_subcortex.smk" +include: "rules/mrtpipelines.smk" diff --git a/dbsc/workflow/rules/mrtpipelines.smk b/dbsc/workflow/rules/mrtpipelines.smk index 7223bc80..66969eed 100644 --- a/dbsc/workflow/rules/mrtpipelines.smk +++ b/dbsc/workflow/rules/mrtpipelines.smk @@ -1,571 +1,523 @@ -from os.path import join +from pathlib import Path +from functools import partial -include: 'zona_bb_subcortex.smk' -include: 'freesurfer.smk' -# DWI Processing -rule nii_to_mif: +include: "zona_bb_subcortex.smk" +include: "freesurfer.smk" + + +# Directories +mrtrix_dir = str(Path(config["output_dir"]) / "mrtrix") + +# Make directory if it doesn't exist +Path(mrtrix_dir).mkdir(parents=True) + +# Parameters +responsemean_flag = config.get("responsemean_dir", None) +shells = config.get("shells", None) +lmax = config.get("lmax", None) + +# BIDS partials +bids_response_out = partial( + bids, + root=mrtrix_dir, + datatype="response", + suffix="response.txt", +) + +bids_dti_out = partial( + bids, + root=mrtrix_dir, + datatype="dti", + model="dti", + **config["subj_wildcards"], +) + +bids_tractography_out = partial( + bids, + root=mrtrix_dir, + datatype="tractography", + **config["subj_wildcards"], +) + +bids_anat_out = partial( + bids, + root=mrtrix_dir, + datatype="anat", + **config["subj_wildcards"], +) + +# Mrtrix3 citation (additional citations are included per rule as necessary): +# Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137 + + +# ------------ MRTRIX PREPROC BEGIN ----------# +rule nii2mif: input: - dwi=bids( - root=join(config["bids_dir"], 'derivatives/prepdwi') - datatype='dwi', - suffix='dwi.nii.gz', - space='T1w', - **config['subj_wildcards'], - ) - bval=bids( - root=join(config["bids_dir"], 'derivatives/prepdwi') - datatype='dwi', - suffix='dwi.bval', - space='T1w', - **config['subj_wildcards'], - ) - bvec=bids( - root=join(config["bids_dir"], 'derivatives/prepdwi') - datatype='dwi', - suffix='dwi.bvec', - space='T1w', - **config['subj_wildcards'], - ) - mask=bids( - root=join(config["bids_dir"], 'derivatives/prepdwi') - datatype='dwi', - suffix='brainmask.nii.gz', - space='T1w', - **config['subj_wildcards'], - ) - params: - threads=workflow.cores, + dwi=config["input_path"]["dwi"], + bval=re.sub(".nii.gz", ".bval", config["input_path"]["dwi"]), + bvec=re.sub(".nii.gz", ".bvec", config["input_path"]["dwi"]), + mask=config["input_path"]["mask"], output: dwi=bids( - root=join(config["output_dir"], 'mrtpipelines'), - datatype='dwi', - suffix='dwi.mif', - space='T1w', - **config['subj_wildcards'] - ) + root=mrtrix_dir, + datatype="dwi", + suffix="dwi.mif", + **config["subj_wildcards"] + ), mask=bids( - root=join(config['output_dir'], 'mrtpipelines'), - datatype='dwi', - suffix='brainmask.mif', - space='T1w', - **config['subj_wildcards'] - ) - group: "subject_1" + root=mrtrix_dir, + datatype="dwi", + suffix="brainmask.mif", + **config["subj_wildcards"] + ), + threads: workflow.cores + group: + "subject_1" container: - config['singularity']['mrtpipelines'] - shell: - 'mrconvert -nthreads {threads} -fslgrad {input.bvecs} {input.bvals} {input.dwi} {output.dwi} ' - 'mrconvert -nthreads {threads} {input.mask} {output.mask}' + config["singularity"]["mrtrix"] + shell: + "mrconvert -nthreads {threads} -fslgrad {input.bvecs} {input.bvals} {input.dwi} {output.dwi} && " + "mrconvert -nthreads {threads} {input.mask} {output.mask}" -rule estimate_response: - """Estimate response functions using Dhollander algorithm """ +rule dwi2response: + """ + Estimate response functions using Dhollander algorithm + + Dhollander, T.; Mito, R.; Raffelt, D. & Connelly, A. Improved white matter response function estimation for 3-tissue constrained spherical deconvolution. Proc Intl Soc Mag Reson Med, 2019, 555 + """ input: - dwi=rules.nii_to_mif.output.dwi, - mask=rules.nii_to_mif.output.mask + dwi=rules.nii2mif.output.dwi, + mask=rules.nii2mif.output.mask, params: - threads=workflow.cores, + shells=f"-shells {shells}" if shells else "", + lmax=f"-lmax {lmax}" if lmax else "", output: - sfwm=bids( - root=join(config['output_dir'], 'mrtpipelines'), - datatype='dwi', - desc='sfwm', - suffix='response.txt', - **config['subj_wildcards'], - ) - gm=bids( - root=join(config['output_dir'], 'mrtpipelines'), - datatype='dwi', - desc='gm' - suffix='response.txt', - **config['subj_wildcards'], - ) - csf=bids( - root=join(config['output_dir'], 'mrtpipelines'), - datatype='dwi', - desc='csf' - suffix='response.txt', - **config['subj_wildcards'], - ) - group: "subject_1" + wm_rf=bids_response_out( + desc="wm", + **config["subj_wildcards"], + ), + gm_rf=bids_response_out( + desc="gm", + **config["subj_wildcards"], + ), + csf_rf=bids_response_out( + desc="csf", + **config["subj_wildcards"], + ), + threads: workflow.cores + group: + "subject_1" container: - config['singularity']['mrtpipelines'] + config["singularity"]["mrtrix"] shell: - 'dwi2response dhollander {input.dwi} {output.sfwm} {output.gm} {output.csf} -nthreads {threads}' + "dwi2response dhollander {input.dwi} {output.wm_rf} {output.gm_rf} {output.csf_rf} -nthreads {threads} -mask {input.mask} {params.shells} {params.lmax} " -# TODO: Add check for pre-computed group response function -rule avg_response: + +rule responsemean: """Compute average response function""" input: - sfwm=expand(rules.estimate_response.output.sfwm, subject=config['subjects']), - gm=expand(rules.estimate_response.output.gm, subject=config['subjects']), - csf=expand(rules.estimate_response.output.csf, subject=config['subjects']) + subject_rf=bids_response_out( + desc="{tissue}", + **config["subj_wildcards"], + ), output: - avg_sfwm=bids( - root=join(config['output_dir'], 'mrtpipelines/avg'), - datatype='dwi', - desc='sfwm', - suffix='response.txt', - ) - avg_gm=bids( - root=join(config['output_dir'], 'mrtpipelines/avg'), - datatype='dwi', - desc='gm', - suffix='response.txt', - ) - avg_csf=bids( - root=join(config['output_dir'], 'mrtpipelines/avg'), - datatype='dwi', - desc='csf', - suffix='response.txt', - ) - group: "group" + avg_rf=bids_response_out( + root=str(Path(mrtrix_dir) / "avg"), + desc="{tissue}", + ), + threads: workflow.cores + group: + "group" container: - config['singularity']['mrtpipelines'] + config["singularity"]["mrtrix"] shell: - 'average_response {input.sfwm} {output.avg_sfwm} {input.gm} {output.avg_gm} {input.csf} {output.csf}' + "responsemean {input.subject_rf} {output.avg_rf} -nthreads {threads}" -rule compute_fod: +rule dwi2fod: + """Jeurissen, B; Tournier, J-D; Dhollander, T; Connelly, A & Sijbers, J. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. NeuroImage, 2014, 103, 411-426""" input: - avg_sfwm=rules.avg_response.output.avg_sfwm, - avg_gm=rules.avg_response.output.avg_gm, - avg_csf=rules.avg_response.output.avg_csf, - dwi=rules.nii_to_mif.output.dwi, - mask=rules.nii_to_mif.output.mask, + dwi=rules.nii2mif.output.dwi, + mask=rules.nii2mif.output.mask, + wm_rf=( + str(Path(config["responsemean_dir"]) / "desc-wm_response.txt") + if responsemean_flag + else expand(rules.responsemean.output.avg_rf, tissue="wm") + ), + gm_rf=( + str(Path(config["responsemean_dir"]) / "desc-gm_response.txt") + if responsemean_flag + else expand(rules.responsemean.output.avg_rf, tissue="gm") + ), + csf_rf=( + str(Path(config["responsemean_dir"]) / "desc-csf_response.txt") + if responsemean_flag + else expand(rules.responsemean.output.csf_avg_rf, tissue="csf") + ), params: - shell=config.get('shells', ''), - lmax=config.get('lmax', ''), - threads=workflow.threads, + shells=f"-shells {shells}" if shells else "", output: - wm_fod=bids( - root=join(config['output_dir'], 'mrtpipelines'), - datatype='fod', - model='csd', - desc='wm', - suffix='fod.mif' - ), - gm_fod=bids( - root=join(config['output_dir'], 'mrtpipelines'), - datatype='fod', - model='csd', - desc='gm', - suffix='fod.mif' - ), - csf_fod=bids( - root=join(config['output_dir'], 'mrtpipelines'), - datatype='fod', - model='csd', - desc='csf', - suffix='fod.mif' - ), - group: "subject_2" + wm_fod=bids_response_out( + model="csd", + desc="wm", + suffix="fod.mif", + **config["subj_wildcards"], + ), + gm_fod=bids_response_out( + model="csd", + desc="gm", + suffix="fod.mif", + **config["subj_wildcards"], + ), + csf_fod=bids_response_out( + model="csd", + desc="csf", + suffix="fod.mif", + **config["subj_wildcards"], + ), + threads: workflow.cores + group: + "subject_2" container: - config['singularity']['mrtpipelines'] + config["singularity"]["mrtrix"] shell: - 'if [ {params.shell} == "" ] && [ {params.lmax} == "" ]; then' - ' dwi2fod -nthreads {params.threads} -mask {input.mask} msmt_csd {input.dwi} {input.avg_sfwm} {output.wm_fod} {input.avg_gm} {output.gm_fod} {input.avg_csf} {output.avg_csf} ' - 'else ' - ' dwi2fod -nthreads {params.threads} -shell {params.shell} -lmax {params.lmax} -mask {input.mask} msmt_csd {input.dwi} {input.avg_sfwm} {output.wm_fod} {input.avg_gm} {output.gm_fod} {input.avg_csf} {output.csf} ' - 'fi' + "dwi2fod -nthreads {threads} {params.shells} -mask {input.mask} msmt_csd {input.dwi} {input.avg_sfwm} {output.wm_fod} {input.avg_gm} {output.gm_fod} {input.avg_csf} {output.csf}" -rule normalise_fod: +rule mtnormalise: + """ + Raffelt, D.; Dhollander, T.; Tournier, J.-D.; Tabbara, R.; Smith, R. E.; Pierre, E. & Connelly, A. Bias Field Correction and Intensity Normalisation for Quantitative Analysis of Apparent Fibre Density. In Proc. ISMRM, 2017, 26, 3541 + + Dhollander, T.; Tabbara, R.; Rosnarho-Tornstrand, J.; Tournier, J.-D.; Raffelt, D. & Connelly, A. Multi-tissue log-domain intensity and inhomogeneity normalisation for quantitative apparent fibre density. In Proc. ISMRM, 2021, 29, 2472 + """ input: - wm_fod=rules.compute_fod.output.wm_fod, - gm_fod=rules.compute_fod.output.gm_fod, - csf_fod=rules.compute_fod.output.csf_fod, - mask=rules.nii_to_mif.output.mask, - params: - threads=workflow.threads, + wm_fod=rules.dwi2fod.output.wm_fod, + gm_fod=rules.dwi2fod.output.gm_fod, + csf_fod=rules.dwi2fod.output.csf_fod, + mask=rules.nii2mif.output.mask, output: - wm_fod=bids( - root=join(config['output_dir'], 'mrtpipelines'), - datatype='fod', - model='csd', - desc='wm', - suffix='fodnorm.mif', - **config['subj_wildcards'], - ), - gm_fod=bids( - root=join(config['output_dir'], 'mrtpipelines'), - datatype='fod', - model='csd', - desc='gm', - suffix='fodnorm.mif', - **config['subj_wildcards'], - ), - csf_fod=bids( - root=join(config['output_dir'], 'mrtpipelines'), - datatype='fod', - model='csd', - desc='csf', - suffix='fodnorm.mif', - **config['subj_wildcards'], - ), - group: "subject_2" + wm_fod=bids_response_out( + model="csd", + desc="wm", + suffix="fodNormalized.mif", + **config["subj_wildcards"], + ), + gm_fod=bids_response_out( + model="csd", + desc="gm", + suffix="fodNormalized.mif", + **config["subj_wildcards"], + ), + csf_fod=bids_response_out( + model="csd", + desc="csf", + suffix="fodNormalized.mif", + **config["subj_wildcards"], + ), + threads: workflow.cores + group: + "subject_2" container: - config['singularity']['mrtpipelines'] + config["singularity"]["mrtrix"] shell: - 'mtnormalise -nthreads {params.threads} -mask {input.mask} {input.wm_fod} {output.wm_fod} {input.gm_fod} {output.gm_fod} {input.csf_fod} {output.csf_fod}' + "mtnormalise -nthreads {threads} -mask {input.mask} {input.wm_fod} {output.wm_fod} {input.gm_fod} {output.gm_fod} {input.csf_fod} {output.csf_fod}" + # DTI (Tensor) Processing -rule dwi_normalise: +rule dwinormalise: input: - dwi=rules.nii_to_mif.output.dwi, - mask=rules.nii_to_mif.output.mask, - params: - threads=workflow.threads, + dwi=rules.nii2mif.output.dwi, + mask=rules.nii2mif.output.mask, output: dwi=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='dwi', - space='T1w', - desc='norm', - suffix='dwi.mif', - **config['subj_wildcards'], - ) + root=mrtrix_dir, + datatype="dwi", + desc="normalized", + suffix="dwi.mif", + **config["subj_wildcards"], + ), + threads: workflow.cores container: - config['singularity']['mrtpipelines'] + config["singularity"]["mrtrix"] + group: + "subject_1" shell: - 'dwinormalise -nthreads {params.threads} {input.dwi} {input.mask} {output.dwi}' - group: participant1 + "dwinormalise individual -nthreads {threads} {input.dwi} {input.mask} {output.dwi}" -rule compute_tensor: +rule dwi2tensor: input: - dwi=rules.dwi_noramlise.dwi, - mask=rules.nii_to_mif.output.mask, - params: - threads=workflow.threads, + dwi=rules.dwinoramlise.dwi, + mask=rules.nii2mif.output.mask, output: - dti=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='dti', - space='T1w', - suffix='dti.mif', - **config['subj_wildcards'], - ), - fa=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='dti', - space='T1w', - model='dti', - suffix='fa.mif', - **config['subj_wildcards'], - ), - ad=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='dti', - space='T1w', - model='dti', - suffix='ad.mif', - **config['subj_wildcards'], - ), - rd=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='dti', - space='T1w', - model='dti', - suffix='rd.mif', - **config['subj_wildcards'], - ), - md=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='dti', - space='T1w', - model='dti', - suffix='fa.mif', - **config['subj_wildcards'], - ) - group: "subject_1" + dti=bids_dti_out(suffix="tensor.mif"), + fa=bids_dti_out(suffix="fa.mif"), + ad=bids_dti_out(suffix="ad.mif"), + rd=bids_dti_out(suffix="rd.mif"), + md=bids_dti_out(suffix="md.mif"), + threads: workflow.cores + group: + "subject_1" container: - config['singularity']['mrtpipelines'] + config["singularity"]["mrtrix"] shell: - 'dwi2tensor -nthreads {params.threads} -mask {input.mask} {input.dwi} {output.dti}' - 'tensor2metric -nthreads {params.threads} -mask {input.mask} {output.dti} -fa {output.fa} -ad {output.ad} -rd {output.rd} -adc {output.md}' + "dwi2tensor -nthreads {threads} -mask {input.mask} {input.dwi} {output.dti} && " + "tensor2metric -nthreads {threads} -mask {input.mask} {output.dti} -fa {output.fa} -ad {output.ad} -rd {output.rd} -adc {output.md}" + +# --------------- MRTRIX PREPROC END --------------# -# Tractography processing -rule gen_tractography: + +# ------------ MRTRIX TRACTOGRAPHY BEGIN ----------# +# TODO (v0.1): CHECK TO MAKE SURE RULES ARE IMPORTED CORRECTLY FROM OTHER SMK FILES + + +rule tckgen: + # Tournier, J.-D.; Calamante, F. & Connelly, A. Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions. Proceedings of the International Society for Magnetic Resonance in Medicine, 2010, 1670 input: - fod=rules.normalise_fod.output.wm_fod, + fod=rules.mtnormalise.output.wm_fod, + mask=rules.nii2mif.output.mask, cortical_ribbon=rules.fs_xfm_to_native.output.ribbon, convex_hull=rules.create_convex_hull.output.convex_hull, subcortical_seg=rules.add_brainstem_new_seg.output.seg, - mask=rules.nii_to_mif.output.mask, params: - threads=workflow.threads, - step=config['step'], - sl=config['sl_count'], + step=config["step"], + sl=config["sl_count"], output: - tck=bids( - root=join(config['out_dir'], mrtpipelines), - datatype='tractography', - space='T1w', - desc='iFOD2', - suffix='tractography.tck', - **config["subj_wildcards"], - ) - group: "subject_2" + tck=bids_tractography_out( + desc="iFOD2", + suffix="tractography.tck", + ), + threads: workflow.cores + group: + "subject_2" container: - config['singularity']['mrtpipelines'] + config["singularity"]["mrtrix"] shell: - 'tckgen -nthreads {params.threads} -algorithm iFOD2 -step {params.step} -select {params.sl_count} -exclude {input.cortical_ribbon} -exclude {input.convex_hull} -include {input.subcortical_seg} -mask {input.mask} -seed_image {input.mask} {input.fod} {output.tck}' + "tckgen -nthreads {threads} -algorithm iFOD2 -step {params.step} -select {params.sl_count} -exclude {input.cortical_ribbon} -exclude {input.convex_hull} -include {input.subcortical_seg} -mask {input.mask} -seed_image {input.mask} {input.fod} {output.tck}" + -rule weight_tractography: +rule tcksift2: + # Smith, R. E.; Tournier, J.-D.; Calamante, F. & Connelly, A. The effects of SIFT on the reproducibility and biological accuracy of the structural connectome. NeuroImage, 2015, 104, 253-265 input: - tck=rules.gen_tractography.output.tck, - fod=rules.normalise_fod.output.wm_fod, - params: - threads=workflow.threads, + tck=rules.tckgen.output.tck, + fod=rules.mtnormalise.output.wm_fod, output: weights=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='iFOD2', - suffix='tckweights.txt', - **config['subj_wildcards'], - ), - mu=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='iFOD2', - suffix='mucoefficient.txt', - **config['subj_wildcards'], - ] - ), - group: "subject_2" + desc="iFOD2", + suffix="tckWeights.txt", + ), + mu=bids(desc="iFOD2", suffix="muCoefficient.txt"), + threads: workflow.cores + group: + "subject_2" container: - config['singularity']['mrtpipelines'] + config["singularity"]["mrtrix"] shell: - 'tcksift2 -nthreads {params.threads} -out_mu {output.mu} {input.tck} {input.fod} {output.weights}' + "tcksift2 -nthreads {threads} -out_mu {output.mu} {input.tck} {input.fod} {output.weights}" -# TODO: ADD OPTION TO OUTPUT TDI MAP -# Connectivity map -# TODO: ADD OPTION TO MULTIPLY BY MU COEFFICIENT -rule connectome_map: +# TODO (v0.2): ADD OPTION TO OUTPUT TDI MAP + + +rule tck2connectome: + """ + Smith, R. E.; Tournier, J.-D.; Calamante, F. & Connelly, A. The effects of SIFT on the reproducibility and biological accuracy of the structural connectome. NeuroImage, 2015, 104, 253-265" + + TODO (v0.2): ADD OPTION TO MULTIPLY BY MU COEFFICIENT + """ input: - weights=rules.weight_tractography.output.weights, - tck=rules.gen_tractography.output.tck, + weights=rules.tcksift2.output.weights, + tck=rules.tckgen.output.tck, subcortical_seg=rules.add_brainstem_new_seg.output.seg, params: - threads=workflow.threads, - radius=config['radial_search'] + radius=config["radial_search"], output: - sl_assignment=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='subcortical', - suffix='nodeassignment.txt', - **config['subj_wildcards'], - ), - node_weights=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='subcortical', - suffix='nodeweights.csv' - **config['subj_wildcards'], - ) - group: "subject_2" + sl_assignment=bids_tractography_out( + desc="subcortical", + suffix="nodeAssignment.txt", + ), + node_weights=bids_tractography_out( + desc="subcortical", + suffix="nodeWeights.csv", + ), + threads: workflow.cores + group: + "subject_2" container: - config['singularity']['mrtpipelines'] + config["singularity"]["mrtrix"] shell: - 'tck2connectome -nthreads {params.threads} -zero_diagonal -stat_edge sum -assignment_radial_search {params.radius} -tck_weights_in {input.weights} -out_assignments {output.sl_assignment} -symmetric {input.tck} {input.subcortical_seg} {output.node_weights} ' - - rule extract_tck: - input: - node_weights=rules.connectome_map.output.node_weights, - sl_assignment=rules.connectome_map.output.sl_assignment, - tck=rules.gen_tractography.output.tck, - params: - threads=workflow.threads, - output: - edge_weight=temp( - bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='subcortical', - suffix='tckweights', - **config['subj_wildcards'], - ) - ), - edge_tck=temp( - bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='subcortical', - suffix='from_', - **config['subj_wildcards'], - ) - ), - group: "subject_2" - container: - config['singularity']['mrtpipelines'], - shell: - 'for i in `seq 2 72`; do ' - ' nodes=$nodes,$i ' - 'done ' - 'connectome2tck -nthreads {params.threads} -nodes $nodes -exclusive -filters_per_edge -tck_weights_in {input.node_weights} -prefix_tck_weights_out {output.edge_weight} {input.tck} {input.sl_assignment} {output.edge_tck} ' - - # NOTE: Pass labelmerge split segs here? - rule create_roi_mask: - # EXPAND OVER NODE1 IN RULE ALL - input: - subcortical_seg=rules.add_brainstem_new_seg.output.seg, - params: - threads=workflow.threads, - output: - roi_mask=temp( - bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='anat', - space='T1w', - desc='{node1}', - suffix='mask.mif' - ) - ), - group: "subject_2" - container: - config['singularity']['mrtpipelines'], - shell: - 'mrcalc -nthreads {params.threads} {input.subcortical_seg} {wildcards.node1} -eq {output.roi_mask}' - - rule filter_tck: - # Node1 should be same as prev rule, need to iterate over node2 - input: - roi1=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='anat', - space='T1w', - desc='{node1}', - suffix='mask.mif' - ), - roi2=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='anat', - space='T1w', - desc='{node2}', - suffix='mask.mif' - ), - # ZI rois (do these need to be separately defined in output?) - lZI=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='anat', - space='T1w', - desc='21', - suffix='mask.mif' - ), - rZI=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='anat', - space='T1w', - desc='22', - suffix='mask.mif' - ), - subcortical_seg=rules.add_brainstem_new_seg.output.seg, - tck=rules.extract_tck.output.edge_tck, - weights=rules.connectome_map.outputs.node_weights - params: - threads=workflow.threads, - output: - filter_mask=temp( - bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='anat', - space='T1w', - desc='exclude', - suffix='mask.mif' - ) - ) - filtered_tck=temp( - bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='from_{node1}-{node2}', - suffix='tractography.tck' - ) + "tck2connectome -nthreads {threads} -zero_diagonal -stat_edge sum -assignment_radial_search {params.radius} -tck_weights_in {input.weights} -out_assignments {output.sl_assignment} -symmetric {input.tck} {input.subcortical_seg} {output.node_weights} " + + +rule connectome2tck: + input: + node_weights=rules.tck2connectome.output.node_weights, + sl_assignment=rules.tck2connectome.output.sl_assignment, + tck=rules.tckgen.output.tck, + params: + nodes=",".join(str(num) for num in range(1, 73)), + output: + edge_weight=temp( + bids_tractography_out( + desc="subcortical", + suffix="tckWeights", ) - filtered_weights=temp( - bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='from_{node1}-{node2}', - suffix='weights.csv' - ) + ), + edge_tck=temp( + bids_tractography_out( + desc="subcortical", + suffix="from", ) - group: "subject_2" - containers: - config['singularity']['mrtpipelines'], - shell: - 'mrcalc -nthreads {params.threads} {input.subcortical_seg} 0 -neq {input.roi1} -sub {input.roi2} -sub {input.lZI} -sub {input.rZI} -sub {output.filter_mask} ' - 'tckedit -nthreads {params.therads} -exclude {output.filter_mask} -tck_weights_in {input.weights} -tck_weights_out {output.filtered_weights} {input.tck} {output.filtered_tck} ' - - rule combine_filtered: - input: - tck=expand(rules.filter_tck.output.filtered_tck, zip, **config['input_zip_lists']['T1w']), - weights=expand(rules.filter_tck.output.filtered_weights, zip, **config['input_zip_lists']['T1w']), - params: - threads=workflow.threads, - output: - combined_tck=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='subcortical', - suffix='tractography.tck' + ), + threads: workflow.cores + group: + "subject_2" + container: + config["singularity"]["mrtrix"] + shell: + "connectome2tck -nthreads {threads} -nodes {params.nodes} -exclusive -filters_per_edge -tck_weights_in {input.node_weights} -prefix_tck_weights_out {output.edge_weight} {input.tck} {input.sl_assignment} {output.edge_tck} " + + +# NOTE: Use labelmerge split segs here? +rule create_roi_mask: + # EXPAND OVER NODE1 IN RULE ALL + input: + subcortical_seg=rules.add_brainstem_new_seg.output.seg, + output: + roi_mask=temp( + bids_anat_out( + desc="{node1}", + suffix="mask.mif", ) - combined_weights=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='subcortical', - suffix='tckweights.txt' + ), + threads: workflow.cores + group: + "subject_2" + container: + config["singularity"]["mrtrix"] + shell: + "mrcalc -nthreads {threads} {input.subcortical_seg} {wildcards.node1} -eq {output.roi_mask}" + + +rule create_exclude_mask: + input: + roi1=bids_anat_out( + desc="{node1}", + suffix="mask.mif", + ), + roi2=bids_anat_out( + desc="{node2}", + suffix="mask.mif", + ), + lZI=bids_anat_out(desc="21", suffix="mask.mif"), + rZI=bids_anat_out(desc="22", suffix="mask.mif"), + subcortical_seg=rules.add_brainstem_new_seg.output.seg, + output: + filter_mask=temp( + bids_anat_out( + desc="exclude{node1}AND{node2}", + suffix="mask.mif", ) - group: "subject_2" - container: - config['singularity']['mrtpipelines'], - shell: - 'tckedit {input.tck} {output.combined_tck} ' - 'cat {input.weights} >> {output.combined_weights} ' - - rule filtered_connectome_map: - input: - weights=rules.combine_filtered.output.combined_weights, - tck=rules.combine_filtered.output.combined_tck, - subcortical_seg=rules.add_brainstem_new_seg.output.seg, - params: - threads=workflow.threads, - radius=config['radial_search'] - output: - sl_assignment=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='subcortical', - suffix='nodeassignment.txt', - **config['subj_wildcards'], + ), + threads: workflow.cores + group: + "subject_2" + container: + config["singularity"]["mrtrix"] + shell: + "mrcalc -nthreads {threads} {input.subcortical_seg} 0 -neq {input.roi1} -sub {input.roi2} -sub {input.lZI} -sub {input.rZI} -sub {output.filter_mask}" + + +rule filter_tck: + input: + filter_mask=rules.create_exclude_mask.outputs.filter_mask, + tck=rules.connectome2tck.output.edge_tck, + weights=rules.tck2connectome.outputs.node_weights, + subcortical_seg=rules.add_brainstem_new_seg.output.seg, + output: + filtered_tck=temp( + bids_tractography_out( + desc="from{node1}to{node2}", + suffix="tractography.tck", ), - node_weights=bids( - root=join(config['out_dir'], 'mrtpipelines'), - datatype='tractography', - space='T1w', - desc='subcortical', - suffix='nodeweights.csv' - **config['subj_wildcards'], - ) - group: "subject_2" - container: - config['singularity']['mrtpipelines'] - shell: - 'tck2connectome -nthreads {params.threads} -zero_diagonal -stat_edge sum -assignment_radial_search {params.radius} -tck_weights_in {input.weights} -out_assignments {output.sl_assignment} -symmetric {input.tck} {input.subcortical_seg} {output.node_weights} -force' \ No newline at end of file + ), + filtered_weights=temp( + bids_tractography_out( + desc="from{node1}to{node2}", + suffix="weights.csv", + ), + ), + threads: workflow.cores + group: + "subject_2" + container: + config["singularity"]["mrtrix"] + shell: + "tckedit -nthreads {threads} -exclude {input.filter_mask} -tck_weights_in {input.weights} -tck_weights_out {output.filtered_weights} {input.tck} {output.filtered_tck}" + + +idxes = np.triu_indices(72, k=1) + + +rule combine_filtered: + input: + tck=expand( + rules.filter_tck.output.filtered_tck, + zip, + node1=list(idxes[0] + 1), + node2=list(idxes[1] + 1), + ), + weights=expand( + rules.filter_tck.output.filtered_weights, + zip, + node1=list(idxes[0] + 1), + node2=list(idxes[1] + 1), + ), + output: + combined_tck=bids_tractography_out( + desc="subcortical", + suffix="tractography.tck", + ), + combined_weights=bids_tractography_out( + desc="subcortical", + suffix="tckWeights.txt", + ), + threads: workflow.cores + group: + "subject_2" + container: + config["singularity"]["mrtrix"] + shell: + "tckedit {input.tck} {output.combined_tck} &&" + "cat {input.weights} >> {output.combined_weights} " + + +rule filtered_tck2connectome: + input: + weights=rules.combine_filtered.output.combined_weights, + tck=rules.combine_filtered.output.combined_tck, + subcortical_seg=rules.add_brainstem_new_seg.output.seg, + params: + radius=config["radial_search"], + output: + sl_assignment=bids_tractography_out( + desc="subcortical", + suffix="nodeAssignment.txt", + ), + node_weights=bids_tractography_out( + desc="subcortical", + suffix="nodeWeights.csv", + ), + threads: workflow.core + group: + "subject_2" + container: + config["singularity"]["mrtrix"] + shell: + "tck2connectome -nthreads {threads} -zero_diagonal -stat_edge sum -assignment_radial_search {params.radius} -tck_weights_in {input.weights} -out_assignments {output.sl_assignment} -symmetric {input.tck} {input.subcortical_seg} {output.node_weights} -force" + + +# ------------ MRTRIX TRACTOGRAPHY END ----------#