General-purpose Snakemake workflow for diffusion-based subcortical parcellation

public 1yr ago 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: diffparc-smk

General-purpose Snakemake workflow for diffusion-based subcortical parcellation

This workflow is customized for the basal forebrain, run on HCP1200 7T diffusion data.

Uses HCP-MMP cortical parcella

Code Snippets

shell: 'fslmaths {input.lh} -max {input.rh} {output.lh_rh} &> {log}'

SnakeMake From line 13 of rules/diffparc.smk

shell:
    'cp {input} {output}'

SnakeMake From line 23 of rules/diffparc.smk

shell:
    'c3d {input} -dilate 1 3x3x3vox -o {output}'

SnakeMake C3D From line 35 of rules/diffparc.smk

shell:
    'antsApplyTransforms -d 3 --interpolation NearestNeighbor -i {input.seed} -o {output} -r {input.ref}  -t {input.invwarp} &> {log}'

SnakeMake From line 56 of rules/diffparc.smk

shell:
    'fslmaths {input.dwi} -bin {output.mask} &&'
    'mri_convert {output.mask} -vs {params.seed_resolution} {params.seed_resolution} {params.seed_resolution} {output.mask_res} -rt nearest &> {log}'

SnakeMake From line 73 of rules/diffparc.smk

shell:
    'reg_resample -flo {input.targets} -res {output.targets_res} -ref {input.mask_res} -NN 0  &> {log}'

SnakeMake From line 90 of rules/diffparc.smk

run:
    f = open(output.target_txt,'w')
    for s in params.target_seg:
        f.write(f'{s}\n')
    f.close()

SnakeMake From line 139 of rules/diffparc.smk

shell:
    'mkdir -p {output.probtrack_dir} && singularity exec -e --nv {params.container} '
    'probtrackx2_gpu --samples={params.bedpost_merged}  --mask={input.mask} --seed={input.seed_res} '
    '--targetmasks={input.target_txt} --seedref={input.seed_res} --nsamples={params.nsamples} '
    '--dir={output.probtrack_dir} {params.probtrack_opts} -V 2  &> {log}'

SnakeMake Singularity Hub From line 167 of rules/diffparc.smk

shell:
    'mkdir -p {output} && ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1 parallel  --jobs {threads} antsApplyTransforms -d 3 --interpolation Linear -i {{1}} -o {{2}}  -r {input.ref} -t {input.warp}  &> {log} :::  {params.in_connmap_3d} :::+ {params.out_connmap_3d}' 

SnakeMake From line 192 of rules/diffparc.smk

script: '../scripts/save_connmap_template_npz.py'

SnakeMake From line 209 of rules/diffparc.smk

script: '../scripts/gather_connmap_group.py'

SnakeMake From line 220 of rules/diffparc.smk

script: '../scripts/spectral_clustering.py'

SnakeMake From line 235 of rules/diffparc.smk

shell:
    '{input.shell_script} {input.seg} {output.seg} {wildcards.k} &> {log}'

SnakeMake From line 248 of rules/diffparc.smk

shell: 
    'mkdir -p {params.out_folder} && tar --extract --file={input.tar} {params.file_in_tar} && '
    'mv {params.file_in_tar} {output.filename}'

SnakeMake From line 26 of rules/hcp_mmp.smk

shell: 'FS_LICENSE={params.license} mris_convert {input} {output} &> {log}'

SnakeMake From line 47 of rules/hcp_mmp.smk

shell: 'FS_LICENSE={params.license} mri_convert {input} {output} &> {log}'

SnakeMake From line 57 of rules/hcp_mmp.smk

shell: 'FS_LICENSE={params.license} mri_info {input.t1} --tkr2scanner > {output.tkr2scanner} 2> {log}'

SnakeMake From line 69 of rules/hcp_mmp.smk

shell: 'wb_command -surface-apply-affine {input.surf} {input.tkr2scanner} {output.surf} &> {log}'

SnakeMake From line 81 of rules/hcp_mmp.smk

shell: 'wb_command -surface-average {output.midthickness} -surf {input.white} -surf {input.pial} &> {log}'

SnakeMake From line 94 of rules/hcp_mmp.smk

shell: 'wb_command -surface-resample {input.surf} {input.current_sphere} {input.new_sphere} {params.method} {output.surf} &> {log}'

SnakeMake From line 111 of rules/hcp_mmp.smk

shell: 
    'wb_command -label-resample {input.label} {input.current_sphere} {input.new_sphere}'
    ' {params.method} {output.label}'
    ' -area-surfs {input.current_surf} {input.new_surf} &> {log}'

SnakeMake From line 130 of rules/hcp_mmp.smk

shell:
    'wb_command -label-to-volume-mapping {input.label} {input.surf} {input.vol_ref} {output.label_vol}'
    ' -ribbon-constrained {input.white_surf} {input.pial_surf}'
    ' -greedy &> {log}'

SnakeMake From line 150 of rules/hcp_mmp.smk

shell:
    'wb_command -label-to-volume-mapping {input.label} {input.surf} {input.vol_ref} {output.label_vol}'
    ' -nearest-vertex {params.nearest_vertex} &> {log}'

SnakeMake From line 172 of rules/hcp_mmp.smk

shell:
    'ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1 parallel  --jobs {threads} antsApplyTransforms -d 3 --interpolation NearestNeighbor -i {{1}} -o {{2}}  -r {input.ref}  -t {input.invwarp} &> {log} :::  {input.cluster_k} :::+ {output.cluster_k}' 

SnakeMake From line 15 of rules/tractmap.smk

shell:
    'fslmaths {input.dwi} -bin {output.mask} &&'
    'mri_convert {output.mask} -vs {params.seed_resolution} {params.seed_resolution} {params.seed_resolution} {output.mask_res} -rt nearest &> {log}'

SnakeMake From line 33 of rules/tractmap.smk

shell:
    'mkdir -p {output.probtrack_dir} && '
    '{params.extract_seed_cmd} && singularity exec -e --nv {params.container} '
    'probtrackx2_gpu --samples={params.bedpost_merged}  --mask={input.mask} --seed={output.probtrack_dir}/in_seed.nii.gz '
    '--seedref={output.probtrack_dir}/in_seed.nii.gz --nsamples={params.nsamples} '
    '--dir={output.probtrack_dir} {params.probtrack_opts} -V 2  &> {log}'

SnakeMake Singularity Hub From line 98 of rules/tractmap.smk

shell:
    'fslmerge -t {output.tractmaps_4d} {input.tractmaps} &> {log}'

SnakeMake From line 118 of rules/tractmap.smk

shell:
    'ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1  antsApplyTransforms -d 3 --interpolation Linear -i {input.tractmap} -o {output.tractmap}  -r {input.ref} -t {input.warp}  &> {log}'

SnakeMake From line 138 of rules/tractmap.smk

shell:
    'fslmerge -t {output.tractmaps_4d} {input.tractmaps} &> {log}'

SnakeMake From line 153 of rules/tractmap.smk

shell:
    'AverageImages 4 {output} 0 {input}'

SnakeMake From line 170 of rules/tractmap.smk

import numpy as np

#load first file to get shape
data = np.load(snakemake.input.connmap_npz[0])
affine = data['affine']
mask = data['mask']
conn_shape = data['conn'].shape
nsubjects = len(snakemake.input.connmap_npz)
conn_group = np.zeros([nsubjects,conn_shape[0],conn_shape[1]])

for i,npz in enumerate(snakemake.input.connmap_npz):
    data = np.load(npz)
    conn_group[i,:,:] = data['conn']

#save conn_group, mask and affine
np.savez(snakemake.output.connmap_group_npz, conn_group=conn_group,mask=mask,affine=affine)

Python numpy From line 1 of scripts/gather_connmap_group.py

import nibabel as nib
import numpy as np
mask_nib = nib.load(snakemake.input.mask)
mask_vol = mask_nib.get_fdata()
mask_indices = mask_vol > 0 
masked = mask_vol[mask_indices]

nvoxels = len(masked)
ntargets = len(snakemake.params.connmap_3d)
conn = np.zeros((nvoxels,ntargets))
for i,conn_file in enumerate(snakemake.params.connmap_3d):
    vol = nib.load(conn_file).get_fdata()
    masked =  vol[mask_indices].T
    conn[:,i] = masked
np.savez(snakemake.output.connmap_npz, conn=conn,mask=mask_vol,affine=mask_nib.affine)

Python numpy nibabel From line 1 of scripts/save_connmap_template_npz.py

import sklearn
import numpy as np
import nibabel as nib

# Define a function for saving niftis 
def save_label_nii (labels,mask,affine,out_nifti):
    labels_vol = np.zeros(mask.shape)
    labels_vol[mask > 0] = labels+1 #add a 1 so label 0 is diff from bgnd
    labels_nib = nib.Nifti1Image(labels_vol,affine)
    nib.save(labels_nib,out_nifti)

data = np.load(snakemake.input.connmap_group_npz)
cluster_range = range(2,snakemake.params.max_k+1)
out_nii_list = snakemake.output

conn_group = data['conn_group']
mask = data['mask']
affine = data['affine']

# Concat subjects
conn_group_m = np.moveaxis(conn_group,0,2)
conn_concat = conn_group_m.reshape([conn_group_m.shape[0],conn_group_m.shape[1]*conn_group_m.shape[2]])

# Run spectral clustering and save output nifti
for i,k in enumerate(cluster_range):
    from sklearn.cluster import SpectralClustering
    clustering = SpectralClustering(n_clusters=k, assign_labels="discretize",random_state=0,affinity='cosine').fit(conn_concat)
    print(f'i={i}, k={k},saving {out_nii_list[i]}')
    save_label_nii(clustering.labels_,mask,affine,out_nii_list[i])