diff --git a/.gitignore b/.gitignore index d79f6bb..984c7ad 100644 --- a/.gitignore +++ b/.gitignore @@ -14,4 +14,7 @@ test_out/ # Misc files .*.swp -report.html \ No newline at end of file +report.html + +# Linters and tools +.ruff_cache \ No newline at end of file diff --git a/snakedwi/resources/qc-app.tar.gz b/snakedwi/resources/qc-app.tar.gz new file mode 100644 index 0000000..4cd5a83 Binary files /dev/null and b/snakedwi/resources/qc-app.tar.gz differ diff --git a/snakedwi/workflow/Snakefile b/snakedwi/workflow/Snakefile index 5e4fe65..6cf0bbc 100644 --- a/snakedwi/workflow/Snakefile +++ b/snakedwi/workflow/Snakefile @@ -27,10 +27,7 @@ include: "rules/diffusion/sdc/synb0.smk" include: "rules/diffusion/eddy.smk" include: "rules/diffusion/reg_dwi_to_t1.smk" include: "rules/diffusion/bedpost.smk" - - -localrules: - check_subj_dwi_metadata, +include: "rules/qc_app.smk" # Currently unused @@ -61,8 +58,8 @@ rule all: input: **get_eddy_quad_all(), **get_bedpost_all(), - mask_qc=expand(rules.qc_b0_brainmask.output.png, zip, **subj_zip_list), - reg_qc=expand(rules.qc_reg_dwi_t1.output.png, zip, **subj_zip_list), + qc=rules.qc.output, + qc_app=rules.unpack_qc_app.output, dtifit=expand( bids( root=root, diff --git a/snakedwi/workflow/lib/check_subj_dwi_metadata.py b/snakedwi/workflow/lib/check_subj_dwi_metadata.py new file mode 100644 index 0000000..edc2b4d --- /dev/null +++ b/snakedwi/workflow/lib/check_subj_dwi_metadata.py @@ -0,0 +1,153 @@ +#!/usr/bin/env python +import json +import sys +from typing import List, Tuple + + +def check_eddy_s2v( + eddy_s2v: bool, + slspec_txt: bool, + has_slice_timing: bool, + json_file: str, +) -> None: + """Check if slice-to-volume correction can be performed""" + if eddy_s2v and not slspec_txt: + if not has_slice_timing: + print( + "Eddy slice-to-volume enabled, but 'SliceTiming' not found in " + f"{json_file}.\n" + "You must do one of the following:\n" + " 1. add the SliceTiming field to your dwi JSON files\n" + " 2. use the '--slspec_txt' option to provide a global slspec " + "file\n" + " 3. do not set the --use_eddy_s2v CLI option to disable " + "slice-to-volume correction" + ) + sys.exit(1) + + +def check_pe_direction(json_dwi: dict, default_pe_dir: str, json_file: str) -> None: + """Check if phase encoding direction is defined""" + if "PhaseEncodingDirection" not in json_dwi: + if not default_pe_dir == "": + print("WARNING: setting default PhaseEncodingDirection") + json_dwi["PhaseEncodingDirection"] = default_pe_dir + else: + if "PhaseEncodingAxis" in json_dwi: + print( + "WARNING: assuming PhaseEncodingDirection from " "PhaseEncodingAxis" + ) + json_dwi["PhaseEncodingDirection"] = json_dwi["PhaseEncodingAxis"] + else: + print( + f"ERROR: PhaseEncodingDirection not found in {json_file}." + "\nYou must add the PhaseEncodingDirection field to your " + "dwi JSON files, or use the " + "--default_phase_encoding_direction CLI option" + ) + sys.exit(1) + + +def check_echo_spacing(json_dwi: dict, default_echo_spacing: str) -> str: + """Check if effective echo spacing is defined""" + if "EffectiveEchoSpacing" in json_dwi: + eff_echo = json_dwi["EffectiveEchoSpacing"] + elif "EstimatedEffectiveEchoSpacing" in json_dwi: + eff_echo = json_dwi["EstimatedEffectiveEchoSpacing"] + else: + print("EffectiveEchoSpacing not defined in JSON, using default value") + eff_echo = default_echo_spacing + + return eff_echo + + +def get_pe_info(json_dwi: dict, pe_axes: List[str], pe_dirs: List[str]) -> None: + """Extract phase encoding information""" + pe_axes.append(json_dwi["PhaseEncodingDirection"][0]) + + if json_dwi["PhaseEncodingDirection"][-1] == "-": + pe_dirs.append("-") + else: + pe_dirs.append("+") + + +def get_pe_indices( + pe_axes: List[str], pe_dirs: List[str], json_files: List[str] +) -> Tuple[List[int], List[str], List[str]]: + """Extract indices to be used in processing""" + # If multiple PE axes use LR/RL if available, otherwise AP otherwise + if len(set(pe_axes)) > 1: + lr_indices = [idx for idx, pe in enumerate(pe_axes) if pe == "i"] + ap_indices = [idx for idx, pe in enumerate(pe_axes) if pe == "j"] + + # If opposite PE directions in LR/RL, otherwise use AP + use_indices = ( + lr_indices + if len(set([pe_dirs[idx] for idx in lr_indices])) == 2 + else ap_indices + ) + else: + use_indices = [idx for idx in range(len(json_files))] + + # Select indices + pe_axes = [pe_axes[idx] for idx in use_indices] + pe_dirs = [pe_dirs[idx] for idx in use_indices] + + return use_indices, pe_axes, pe_dirs + + +def check_subj_dwi_metadata( + json_files: List[str], + smk_config: dict, +) -> None: + bool_map = {True: "yes", False: "no"} + pe_axes, pe_dirs = [], [] + eddy_s2v = smk_config["use_eddy_s2v"] + slspec_txt = smk_config["slspec_txt"] + + for json_file in json_files: + with open(json_file, "r") as fname: + json_dwi = json.load(fname) + + has_slice_timing = True if "SliceTiming" in json_dwi else False + + # Slice-to-volume + check_eddy_s2v( + eddy_s2v=eddy_s2v, + slspec_txt=slspec_txt, + has_slice_timing=has_slice_timing, + json_file=json_file, + ) + # Phase encoding direction + check_pe_direction( + json_dwi=json_dwi, + default_pe_dir=smk_config["default_phase_encoding_direction"], + json_file=json_file, + ) + # Effective echo spacing (not used) + eff_echo = check_echo_spacing( # noqa: F841 + json_dwi=json_dwi, + default_echo_spacing=smk_config["default_effective_echo_spacing"], + ) + # Save PE information + get_pe_info(json_dwi=json_dwi, pe_axes=pe_axes, pe_dirs=pe_dirs) + + # Get PE indices + use_indices, pe_axes, pe_dirs = get_pe_indices( + pe_axes=pe_axes, pe_dirs=pe_dirs, json_files=json_files + ) + + # Indices + settings = { + "indices": use_indices, + "eddys2v": bool_map[eddy_s2v], + "PEaxis": pe_axes[0], + } + if smk_config["sdc_method"] == "optimal": + settings["sdc"] = ( + "topup" if len(set(pe_dirs)) >= 2 else smk_config["sdc_method_alternate"] + ) + else: + settings["sdc"] = smk_config["sdc_method"] + + return settings diff --git a/snakedwi/workflow/rules/diffusion/method_grabber.smk b/snakedwi/workflow/rules/diffusion/method_grabber.smk index 59d2632..738c671 100644 --- a/snakedwi/workflow/rules/diffusion/method_grabber.smk +++ b/snakedwi/workflow/rules/diffusion/method_grabber.smk @@ -53,10 +53,7 @@ def get_dwi_ref(wildcards): **subj_wildcards, ) else: - checkpoint_output = checkpoints.check_subj_dwi_metadata.get(**wildcards).output[ - 0 - ] - ([method],) = glob_wildcards(os.path.join(checkpoint_output, "sdc-{method}")) + method = get_sdc_method(wildcards) return sdc_methods[method] diff --git a/snakedwi/workflow/rules/diffusion/prepdwi.smk b/snakedwi/workflow/rules/diffusion/prepdwi.smk index 9314326..75f7d93 100644 --- a/snakedwi/workflow/rules/diffusion/prepdwi.smk +++ b/snakedwi/workflow/rules/diffusion/prepdwi.smk @@ -458,3 +458,19 @@ rule qc_b0_brainmask: config["singularity"]["python"] script: "../../scripts/qc/vis_qc_dseg.py" + + +rule compile_qc_b0_brainmask_manifest: + input: + mask_qc=expand(rules.qc_b0_brainmask.output.png, zip, **subj_zip_list), + output: + os.path.join(qc, "data", "mask.json"), + run: + with open(output[0], "w") as f: + json.dump( + { + "title": "DWI Brainmask", + "images": sorted(input["mask_qc"]), + }, + f, + ) diff --git a/snakedwi/workflow/rules/diffusion/reg_dwi_to_t1.smk b/snakedwi/workflow/rules/diffusion/reg_dwi_to_t1.smk index 68598e6..bdd3961 100644 --- a/snakedwi/workflow/rules/diffusion/reg_dwi_to_t1.smk +++ b/snakedwi/workflow/rules/diffusion/reg_dwi_to_t1.smk @@ -175,6 +175,22 @@ rule qc_reg_dwi_t1: "../../scripts/qc/vis_regqc.py" +rule compile_qc_reg_dwi_t1_manifest: + input: + reg_qc=expand(rules.qc_reg_dwi_t1.output.png, zip, **subj_zip_list), + output: + os.path.join(qc, "data", "reg.json"), + run: + with open(output[0], "w") as f: + json.dump( + { + "title": "T1w Registration", + "images": sorted(input["reg_qc"]), + }, + f, + ) + + rule convert_xfm_ras2itk: input: xfm_ras=rules.reg_dwi_to_t1.output.xfm_ras, diff --git a/snakedwi/workflow/rules/diffusion/sdc/synthsr.smk b/snakedwi/workflow/rules/diffusion/sdc/synthsr.smk index 854c971..69b37dd 100644 --- a/snakedwi/workflow/rules/diffusion/sdc/synthsr.smk +++ b/snakedwi/workflow/rules/diffusion/sdc/synthsr.smk @@ -97,7 +97,6 @@ def get_restrict_deformation(wildcards, input, output): rule reg_b0_to_t1_synthsr: input: - rules.check_subj_dwi_metadata.output, t1synth=rules.rigid_reg_t1_to_b0_synthsr.output.warped_subj, b0synth=rules.reslice_synthSR_b0.output.synthsr, params: diff --git a/snakedwi/workflow/rules/diffusion/sdc/topup.smk b/snakedwi/workflow/rules/diffusion/sdc/topup.smk index 3c4a4fd..46ff720 100644 --- a/snakedwi/workflow/rules/diffusion/sdc/topup.smk +++ b/snakedwi/workflow/rules/diffusion/sdc/topup.smk @@ -80,9 +80,6 @@ rule apply_topup_jac: phenc_concat=rules.concat_phase_encode_txt.output.phenc_concat, topup_fieldcoef=rules.run_topup.output.topup_fieldcoef, topup_movpar=rules.run_topup.output.topup_movpar, - workflowopts=bids( - root=root, datatype="dwi", suffix="workflowopts", **subj_wildcards - ), params: inindex=get_applytopup_inindex, topup_prefix=bids( diff --git a/snakedwi/workflow/rules/gradcorrect.smk b/snakedwi/workflow/rules/gradcorrect.smk index 0556cb8..d8a8ea6 100644 --- a/snakedwi/workflow/rules/gradcorrect.smk +++ b/snakedwi/workflow/rules/gradcorrect.smk @@ -43,10 +43,7 @@ sdc_methods = { def get_dwi_ref_for_gradcorrect(wildcards): if config["gradcorrect_coeffs"]: - checkpoint_output = checkpoints.check_subj_dwi_metadata.get(**wildcards).output[ - 0 - ] - ([method],) = glob_wildcards(os.path.join(checkpoint_output, "sdc-{method}")) + method = get_sdc_method(wildcards) return sdc_methods[method] else: diff --git a/snakedwi/workflow/rules/qc_app.smk b/snakedwi/workflow/rules/qc_app.smk new file mode 100644 index 0000000..c271a0b --- /dev/null +++ b/snakedwi/workflow/rules/qc_app.smk @@ -0,0 +1,38 @@ +rule qc: + input: + mask_qc=rules.compile_qc_b0_brainmask_manifest.output[0], + reg_qc=rules.compile_qc_reg_dwi_t1_manifest.output[0], + output: + os.path.join(qc, "data.json"), + run: + with open(output[0], "w") as f: + json.dump( + { + "mask": json.loads(Path(input["mask_qc"]).read_text()), + "reg": json.loads(Path(input["reg_qc"]).read_text()), + }, + f, + ) + + +_qc_app = os.path.join(workflow.basedir, "..", "resources", "qc-app.tar.gz") + + +def _get_tar_contents(file): + try: + return [ + p + for p in sp.check_output(["tar", "-tf", _qc_app]).decode().splitlines() + if p[-1] != "/" + ] + except sp.CalledProcessError as err: + raise Exception("Unable to find qc-app.tar.gz...") from err + + +rule unpack_qc_app: + input: + os.path.join(workflow.basedir, "..", "resources", "qc-app.tar.gz"), + output: + _get_tar_contents(_qc_app), + shell: + "tar -xvzf {input}" diff --git a/snakedwi/workflow/rules/setup.smk b/snakedwi/workflow/rules/setup.smk index be626af..6bc5146 100644 --- a/snakedwi/workflow/rules/setup.smk +++ b/snakedwi/workflow/rules/setup.smk @@ -1,3 +1,5 @@ +from functools import lru_cache +import json from pathlib import Path from snakebids import ( bids, @@ -5,8 +7,10 @@ from snakebids import ( filter_list, get_wildcard_constraints, ) +import subprocess as sp # from snakeboost import PipEnv +from lib.check_subj_dwi_metadata import check_subj_dwi_metadata # writes inputs_config.yml and updates config dict @@ -39,6 +43,7 @@ input_path = inputs.input_path root = os.path.join(config["root"], "snakedwi") work = os.path.join(config["root"], "work") +qc = os.path.join(config["root"], "qc") # setup pipenvs - all my python rules use the script: directive, so will be some work to use snakeboost for this.. # dwi_env = PipEnv( # root=Path('work'), diff --git a/snakedwi/workflow/rules/workflowopts.smk b/snakedwi/workflow/rules/workflowopts.smk index 176eaed..1af6e6f 100644 --- a/snakedwi/workflow/rules/workflowopts.smk +++ b/snakedwi/workflow/rules/workflowopts.smk @@ -1,4 +1,16 @@ -checkpoint check_subj_dwi_metadata: +@lru_cache(None) +def _get_settings(**wildcards): + return check_subj_dwi_metadata( + json_files=expand( + re.sub(".nii.gz", ".json", input_path["dwi"]), + zip, + **filter_list(input_zip_lists["dwi"], wildcards) + ), + smk_config=config, + ) + + +rule check_subj_dwi_metadata: input: dwi_jsons=lambda wildcards: expand( re.sub(".nii.gz", ".json", input_path["dwi"]), @@ -13,14 +25,6 @@ checkpoint check_subj_dwi_metadata: ), index_col_name="subj", output: - workflowopts=directory( - bids( - root=root, - datatype="dwi", - suffix="workflowopts", - **subj_wildcards, - ) - ), metadata=bids( root=root, datatype="dwi", @@ -59,62 +63,36 @@ rule create_missing_subj_tsv: def get_dwi_indices(all_dwi, wildcards): - checkpoint_output = checkpoints.check_subj_dwi_metadata.get(**wildcards).output[0] - ([indices_str],) = glob_wildcards( - os.path.join(checkpoint_output, "indices-{indices}") - ) - indices = indices_str.split(",") - return [all_dwi[int(i)] for i in indices] + return [all_dwi[i] for i in _get_settings(**wildcards)["indices"]] def get_dwi_indices_only(wildcards): - checkpoint_output = checkpoints.check_subj_dwi_metadata.get(**wildcards).output[0] - ([indices_str],) = glob_wildcards( - os.path.join(checkpoint_output, "indices-{indices}") - ) - indices = indices_str.split(",") - return [int(i) for i in indices] + return _get_settings(**wildcards)["indices"] def get_sdc_method(wildcards): - checkpoint_output = checkpoints.check_subj_dwi_metadata.get(**wildcards).output[0] - # this gets the sdc method for this subject(session) - ([method],) = glob_wildcards(os.path.join(checkpoint_output, "sdc-{method}")) - return method + return _get_settings(**wildcards)["sdc"] def get_dwi_num_scans(wildcards): # this gets the number of DWI scans for this subject(session) - checkpoint_output = checkpoints.check_subj_dwi_metadata.get(**wildcards).output[0] - ([indices_str],) = glob_wildcards( - os.path.join(checkpoint_output, "indices-{indices}") - ) - indices = indices_str.split(",") - return len(indices) + return len(_get_settings(**wildcards)["indices"]) def get_pe_axis(wildcards): - checkpoint_output = checkpoints.check_subj_dwi_metadata.get(**wildcards).output[0] - ([peaxis],) = glob_wildcards(os.path.join(checkpoint_output, "PEaxis-{axis}")) - return peaxis + return _get_settings(**wildcards)["PEaxis"] def get_enable_s2v(wildcards): - checkpoint_output = checkpoints.check_subj_dwi_metadata.get(**wildcards).output[0] - ([s2v_is_enabled],) = glob_wildcards( - os.path.join(checkpoint_output, "eddys2v-{isenabled}") - ) - return s2v_is_enabled + return _get_settings(**wildcards)["eddys2v"] def get_index_of_dwi_scan(wildcards): """given wildcards into a specific dwi acquisition, this returns the 0-based index of that scan in the list of dwi scans in use for this acquisition""" - checkpoint_output = checkpoints.check_subj_dwi_metadata.get(**wildcards).output[0] - ([indices_str],) = glob_wildcards( - os.path.join(checkpoint_output, "indices-{indices}") + dwi_indices = get_dwi_indices_only( + {wcard: wildcards[wcard] for wcard in inputs.subj_wildcards} ) - dwi_indices = [int(i) for i in indices_str.split(",")] # first filter the ziplist to get the subject(+session) subj_filter = {"subject": wildcards.subject} diff --git a/snakedwi/workflow/scripts/metadata/check_subj_dwi_metadata.py b/snakedwi/workflow/scripts/metadata/check_subj_dwi_metadata.py index d9219b7..cc63697 100644 --- a/snakedwi/workflow/scripts/metadata/check_subj_dwi_metadata.py +++ b/snakedwi/workflow/scripts/metadata/check_subj_dwi_metadata.py @@ -4,7 +4,6 @@ from typing import List, Tuple import pandas as pd -from snakemake.shell import shell def check_eddy_s2v( @@ -99,26 +98,9 @@ def get_pe_indices( return use_indices, pe_axes, pe_dirs -def set_sdc_method( - workflowopts: str, - pe_dirs: List[str], - smk_config: dict, -) -> None: - # if using "optimal" option vs user choice (e.g. topup, synthsr) - if smk_config["sdc_method"] == "optimal": - if len(set(pe_dirs)) >= 2: # Opp PE - shell(f"touch {workflowopts}/sdc-topup") - else: # Single PE - shell(f"touch {workflowopts}/sdc-{smk_config['sdc_method_alternate']}") - - else: - shell(f"touch {workflowopts}/sdc-{smk_config['sdc_method']}") - - def check_subj_dwi_metadata( json_files: List[str], index_col: Tuple[str, str], - workflowopts: str, metadata: dict, smk_config: dict, ) -> None: @@ -160,16 +142,6 @@ def check_subj_dwi_metadata( ) # Write workflow opts (as workflow trigger): - shell(f"mkdir -p {workflowopts}") - # Indices - write_indices = ",".join([f"{idx}" for idx in use_indices]) - shell(f"touch {workflowopts}/indices-{write_indices}") - # Set SDC method - set_sdc_method(workflowopts=workflowopts, pe_dirs=pe_dirs, smk_config=smk_config) - # Slice-to-volume correction - shell(f"touch {workflowopts}/eddys2v-{bool_map[eddy_s2v]}") - # PE Axis - shell(f"touch {workflowopts}/PEaxis-{pe_axes[0]}") # Write metadata dict metadata_dict = { @@ -190,7 +162,6 @@ def check_subj_dwi_metadata( snakemake.params.index_col_value, # noqa: F821 snakemake.params.index_col_name, # noqa: F821 ), - workflowopts=snakemake.output.workflowopts, # noqa: F821 metadata=snakemake.output.metadata, # noqa: F821 smk_config=snakemake.config, # noqa: F821 )