diff --git a/README.md b/README.md
index e5d5b8d..018def7 100644
--- a/README.md
+++ b/README.md
@@ -18,6 +18,7 @@ Takes TIF images (tiled or prestitched) and outputs a validated BIDS Microscopy
- Python >= 3.11
- Lightsheet data:
- Raw Ultramicroscope Blaze OME TIFF files (include `blaze` in the acquisition tag)
+ - can be 2D or 3D TIFF files
- Prestitched TIFF files (include `prestitched` in the acquisition tag)
@@ -63,10 +64,8 @@ or for snakemake<8.0, use:
snakemake -c all --use-singularity
```
-Note: if you run the workflow on a system with large memory, you will need to set the heap size for the stitching and fusion rules. This can be done with e.g.: `--set-resources bigstitcher_spark_stitching:mem_mb=60000 bigstitcher_spark_fusion:mem_mb=100000`
+Note: if you run the workflow on a system with large memory, you will need to set the heap size for the stitching and fusion rules. This can be done with e.g.: `--set-resources bigstitcher_stitching:mem_mb=60000 bigstitcher_fusion:mem_mb=100000`
7. If you want to run the workflow using a batch job submission server, please see the executor plugins here: https://snakemake.github.io/snakemake-plugin-catalog/
-
-Alternate usage of this workflow (making use of conda) is described in the [Snakemake Workflow Catalog](https://snakemake.github.io/snakemake-workflow-catalog?repo=khanlab/SPIMprep).
diff --git a/config/config.yml b/config/config.yml
index edc16ba..9858eb7 100644
--- a/config/config.yml
+++ b/config/config.yml
@@ -1,6 +1,5 @@
datasets: 'config/datasets.tsv'
-
root: 'bids' # can use a s3:// or gcs:// prefix to write output to cloud storage
work: 'work'
@@ -13,6 +12,7 @@ cores_per_rule: 32
#import wildcards: tilex, tiley, channel, zslice (and prefix - unused)
import_blaze:
raw_tif_pattern: "{prefix}_Blaze[{tilex} x {tiley}]_C{channel}_xyz-Table Z{zslice}.ome.tif"
+ raw_tif_pattern_zstack: "{prefix}_Blaze[{tilex} x {tiley}]_C{channel}.ome.tif"
intensity_rescaling: 0.5 #raw images seem to be at the upper end of uint16 (over-saturated) -- causes wrapping issues when adjusting with flatfield correction etc. this rescales the raw data as it imports it..
import_prestitched:
@@ -36,12 +36,10 @@ bigstitcher:
downsample_in_x: 4
downsample_in_y: 4
downsample_in_z: 1
- method: "phase_corr" #unused
- methods: #unused
+ methods: #unused, only for reference
phase_corr: "Phase Correlation"
optical_flow: "Lucas-Kanade"
filter_pairwise_shifts:
- enabled: 1 #unused
min_r: 0.7
max_shift_total: 50
global_optimization:
@@ -64,7 +62,7 @@ bigstitcher:
block_size_factor_z: 32
ome_zarr:
- desc: sparkstitchedflatcorr
+ desc: stitchedflatcorr
max_downsampling_layers: 5 # e.g. 4 levels: { 0: orig, 1: ds2, 2: ds4, 3: ds8, 4: ds16}
rechunk_size: #z, y, x
- 1
@@ -155,5 +153,5 @@ report:
containers:
- spimprep: 'docker://khanlab/spimprep-deps:main'
+ spimprep: 'docker://khanlab/spimprep-deps:v0.1.0'
diff --git a/resources/qc/ff_html_temp.html b/resources/qc/ff_html_temp.html
index 658b26d..ba6ee7a 100644
--- a/resources/qc/ff_html_temp.html
+++ b/resources/qc/ff_html_temp.html
@@ -43,13 +43,13 @@ Chunk - {{ chunk_num }} Channel - {{ loop.index0 }}
{%- for image in channel %}
- 
- Corrected
+ 
+ Uncorrected
Slice-{{ image.slice }}
|
-
- Uncorrected
+
+ Corrected
Slice-{{ image.slice }}
|
{%- endfor %}
@@ -117,4 +117,4 @@ Uncorrected
}
})
-
\ No newline at end of file
+
diff --git a/workflow/Snakefile b/workflow/Snakefile
index cd19ab3..c997f26 100644
--- a/workflow/Snakefile
+++ b/workflow/Snakefile
@@ -33,7 +33,7 @@ rule all:
input:
get_all_targets(),
get_bids_toplevel_targets(),
- get_qc_targets(), #need to skip this if using prestitched
+# get_qc_targets(), #need to skip this if using prestitched
localrule: True
diff --git a/workflow/lib/dask_image.py b/workflow/lib/dask_image.py
new file mode 100644
index 0000000..dbc06c3
--- /dev/null
+++ b/workflow/lib/dask_image.py
@@ -0,0 +1,45 @@
+from __future__ import annotations
+
+import os
+from glob import glob
+
+try:
+ import tifffile
+except (AttributeError, ImportError):
+ pass
+
+from dask.array.core import Array
+from dask.base import tokenize
+
+
+def add_leading_dimension(x):
+ return x[None, ...]
+
+def imread(fn, page):
+ return tifffile.imread(fn,key=page)
+
+
+def imread_pages(filename, preprocess=None):
+
+
+ tif = tifffile.TiffFile(filename)
+ pages = [i for i in range(len(tif.pages))]
+ name = "imread-%s" % tokenize(filename, os.path.getmtime(filename))
+
+ sample = tif.pages[0].asarray()
+
+ if preprocess:
+ sample = preprocess(sample)
+
+ keys = [(name, i) + (0,) * len(sample.shape) for i in pages]
+ if preprocess:
+ values = [
+ (add_leading_dimension, (preprocess, (imread, filename,i))) for i in pages
+ ]
+ else:
+ values = [(add_leading_dimension, (imread, filename,i)) for i in pages]
+ dsk = dict(zip(keys, values))
+
+ chunks = ((1,) * len(pages),) + tuple((d,) for d in sample.shape)
+
+ return Array(dsk, name, chunks, sample.dtype)
diff --git a/workflow/macros/AutostitchMacro.ijm b/workflow/macros/AutostitchMacro.ijm
deleted file mode 100644
index e9df4df..0000000
--- a/workflow/macros/AutostitchMacro.ijm
+++ /dev/null
@@ -1,56 +0,0 @@
-args = getArgument()
-args = split(args, " ");
-
-dataset_xml = args[0]
-pairwise_method = args[1]
-ds_x = args[2]
-ds_y = args[3]
-ds_z = args[4]
-do_filter = args[5]
-min_r = args[6]
-do_global = args[7]
-global_strategy = args[8]
-
-run("Calculate pairwise shifts ...",
- "select=" + dataset_xml +
- " process_angle=[All angles] " +
- " process_channel=[All channels] " +
- " process_illumination=[All illuminations] " +
- " process_tile=[All tiles] " +
- " process_timepoint=[All Timepoints] " +
- " method=["+ pairwise_method +"] " +
- " channels=[Average Channels] " +
- " downsample_in_x=" + ds_x +
- " downsample_in_y=" + ds_y +
- " downsample_in_z=" + ds_z );
-
-
-if ( do_filter == 1 ){
-
-run("Filter pairwise shifts ...",
- "select=" + dataset_xml +
- " min_r=" + min_r +
- " max_r=1 " +
- " max_shift_in_x=0 " +
- " max_shift_in_y=0 " +
- " max_shift_in_z=0 " +
- " max_displacement=0");
-}
-
-if ( do_global == 1 ){
-
-run("Optimize globally and apply shifts ...",
- "select=" + dataset_xml +
- " process_angle=[All angles] " +
- " process_channel=[All channels] " +
- " process_illumination=[All illuminations] " +
- " process_tile=[All tiles] " +
- " process_timepoint=[All Timepoints] " +
- " relative=2.500 " +
- " absolute=3.500 " +
- " global_optimization_strategy=["+global_strategy+"] " +
- " fix_group_0-0,");
-}
-
-// quit after we are finished
-eval("script", "System.exit(0);");
diff --git a/workflow/macros/FuseImageMacroZarr.ijm b/workflow/macros/FuseImageMacroZarr.ijm
deleted file mode 100644
index 9d49b06..0000000
--- a/workflow/macros/FuseImageMacroZarr.ijm
+++ /dev/null
@@ -1,47 +0,0 @@
-args = getArgument()
-args = split(args, " ");
-
-in_dataset_xml = args[0]
-downsampling = args[1]
-channel = args[2]
-output_zarr = args[3]
-block_size_x = args[4]
-block_size_y = args[5]
-block_size_z = args[6]
-block_size_factor_x = args[7]
-block_size_factor_y = args[8]
-block_size_factor_z = args[9]
-
-run("Fuse dataset ...",
- "select=" + in_dataset_xml +
- " process_angle=[All angles] " +
- " process_channel=[Single channel (Select from List)] " +
- " processing_channel=[channel " + channel + "]" +
- " process_illumination=[All illuminations] " +
- " process_tile=[All tiles] " +
- " process_timepoint=[All Timepoints] " +
- " bounding_box=[Currently Selected Views] " +
- " preserve_original " +
- " downsampling=" + downsampling +
- " interpolation=[Linear Interpolation] " +
- " pixel_type=[16-bit unsigned integer] " +
- " interest_points_for_non_rigid=[-= Disable Non-Rigid =-] " +
- " blend produce=[Each timepoint & channel] " +
- " fused_image=[ZARR/N5/HDF5 export using N5-API] " +
- " define_input=[Auto-load from input data (values shown below)] " +
- " export=ZARR create_0 " +
- " zarr_dataset_path=" + output_zarr +
- " zarr_base_dataset=/ " +
- " zarr_dataset_extension=/s0 " +
- " show_advanced_block_size_options " +
- " block_size_x=" + block_size_x +
- " block_size_y=" + block_size_y +
- " block_size_z=" + block_size_z +
- " block_size_factor_x=" + block_size_factor_x +
- " block_size_factor_y=" + block_size_factor_y +
- " block_size_factor_z=" + block_size_factor_z +
- " subsampling_factors=[{ {1,1,1}}]");
-
-
-// quit after we are finished
-eval("script", "System.exit(0);");
diff --git a/workflow/rules/bigstitcher.smk b/workflow/rules/bigstitcher.smk
index 4256e3b..24a4a3e 100644
--- a/workflow/rules/bigstitcher.smk
+++ b/workflow/rules/bigstitcher.smk
@@ -94,74 +94,7 @@ rule zarr_to_bdv:
"../scripts/zarr_to_n5_bdv.py"
-rule bigstitcher:
- input:
- dataset_n5=rules.zarr_to_bdv.output.bdv_n5,
- dataset_xml=rules.zarr_to_bdv.output.bdv_xml,
- ijm=Path(workflow.basedir) / "macros" / "AutostitchMacro.ijm",
- params:
- fiji_launcher_cmd=get_fiji_launcher_cmd,
- macro_args=get_macro_args_bigstitcher,
- rm_old_xml=lambda wildcards, output: f"rm -f {output.dataset_xml}~?",
- output:
- launcher=temp(
- bids(
- root=work,
- subject="{subject}",
- datatype="micr",
- sample="{sample}",
- acq="{acq}",
- desc="{desc}",
- suffix="bigstitcherfiji.sh",
- )
- ),
- dataset_xml=temp(
- bids(
- root=work,
- subject="{subject}",
- datatype="micr",
- sample="{sample}",
- acq="{acq}",
- desc="{desc}",
- suffix="bigstitcherfiji.xml",
- )
- ),
- benchmark:
- bids(
- root="benchmarks",
- datatype="bigstitcherfiji",
- subject="{subject}",
- sample="{sample}",
- acq="{acq}",
- desc="{desc}",
- suffix="benchmark.tsv",
- )
- log:
- bids(
- root="logs",
- datatype="bigstitcherfiji",
- subject="{subject}",
- sample="{sample}",
- acq="{acq}",
- desc="{desc}",
- suffix="log.txt",
- ),
- container:
- config["containers"]["spimprep"]
- resources:
- runtime=30,
- mem_mb=10000,
- threads: config["cores_per_rule"]
- group:
- "preproc"
- shell:
- "cp {input.dataset_xml} {output.dataset_xml} && "
- " {params.fiji_launcher_cmd} && "
- " echo ' -macro {input.ijm} \"{params.macro_args}\"' >> {output.launcher} "
- " && {output.launcher} |& tee {log} && {params.rm_old_xml}"
-
-
-rule bigstitcher_spark_stitching:
+rule bigstitcher_stitching:
input:
dataset_n5=rules.zarr_to_bdv.output.bdv_n5,
dataset_xml=rules.zarr_to_bdv.output.bdv_xml,
@@ -230,10 +163,10 @@ rule bigstitcher_spark_stitching:
"{params.rm_old_xml}"
-rule bigstitcher_spark_solver:
+rule bigstitcher_solver:
input:
dataset_n5=rules.zarr_to_bdv.output.bdv_n5,
- dataset_xml=rules.bigstitcher_spark_stitching.output.dataset_xml,
+ dataset_xml=rules.bigstitcher_stitching.output.dataset_xml,
params:
downsampling="--downsampling={dsx},{dsy},{dsz}".format(
dsx=config["bigstitcher"]["calc_pairwise_shifts"]["downsample_in_x"],
@@ -295,7 +228,6 @@ rule bigstitcher_spark_solver:
"{params.rm_old_xml}"
#lambda 0.1 is default (can expose this if needed)
-
rule bigstitcher_fusion:
input:
dataset_n5=bids(
@@ -320,98 +252,6 @@ rule bigstitcher_fusion:
else "stitching"
),
),
- ijm=Path(workflow.basedir) / "macros" / "FuseImageMacroZarr.ijm",
- params:
- fiji_launcher_cmd=get_fiji_launcher_cmd,
- macro_args=get_macro_args_zarr_fusion,
- output:
- launcher=temp(
- bids(
- root=work,
- subject="{subject}",
- datatype="micr",
- sample="{sample}",
- acq="{acq}",
- desc="stitched{desc}",
- stain="{stain}",
- suffix="fuseimagen5.sh",
- )
- ),
- zarr=temp(
- directory(
- bids(
- root=work,
- subject="{subject}",
- datatype="micr",
- sample="{sample}",
- acq="{acq}",
- desc="stitched{desc}",
- stain="{stain}",
- suffix="SPIM.zarr",
- )
- )
- ),
- benchmark:
- bids(
- root="benchmarks",
- datatype="fuse_dataset",
- subject="{subject}",
- sample="{sample}",
- acq="{acq}",
- desc="stitched{desc}",
- stain="{stain}",
- suffix="benchmark.tsv",
- )
- log:
- bids(
- root="logs",
- datatype="fuse_dataset",
- subject="{subject}",
- sample="{sample}",
- acq="{acq}",
- desc="stitched{desc}",
- stain="{stain}",
- suffix="log.txt",
- ),
- container:
- config["containers"]["spimprep"]
- resources:
- runtime=30,
- mem_mb=40000,
- threads: config["cores_per_rule"]
- group:
- "preproc"
- shell:
- " {params.fiji_launcher_cmd} && "
- " echo ' -macro {input.ijm} \"{params.macro_args}\"' >> {output.launcher} "
- " && {output.launcher} |& tee {log}"
-
-
-rule bigstitcher_spark_fusion:
- input:
- dataset_n5=bids(
- root=work,
- subject="{subject}",
- datatype="micr",
- sample="{sample}",
- acq="{acq}",
- desc="{desc}",
- suffix="bdv.n5",
- ),
- dataset_xml=bids(
- root=work,
- subject="{subject}",
- datatype="micr",
- sample="{sample}",
- acq="{acq}",
- desc="{desc}",
- suffix="bigstitcher{}.xml".format(
- "solver"
- if config["bigstitcher"]["global_optimization"]["enabled"]
- else "stitching"
- ),
- ),
- ijm=Path(workflow.basedir) / "macros" / "FuseImageMacroZarr.ijm",
params:
channel=lambda wildcards: "--channelId={channel}".format(
channel=get_stains(wildcards).index(wildcards.stain)
@@ -438,7 +278,7 @@ rule bigstitcher_spark_fusion:
datatype="micr",
sample="{sample}",
acq="{acq}",
- desc="sparkstitched{desc}",
+ desc="stitched{desc}",
stain="{stain}",
suffix="SPIM.zarr",
)
diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk
index 88d40e7..80a773a 100644
--- a/workflow/rules/common.smk
+++ b/workflow/rules/common.smk
@@ -257,52 +257,6 @@ def get_stains(wildcards):
return df.iloc[0][stain_columns].dropna().tolist()
-# bigstitcher
-def get_fiji_launcher_cmd(wildcards, output, threads, resources):
- launcher_opts_find = "-Xincgc"
- launcher_opts_replace = f"-XX:+UseG1GC -verbose:gc -XX:+PrintGCDateStamps -XX:ActiveProcessorCount={threads}"
- pipe_cmds = []
- pipe_cmds.append("ImageJ-linux64 --dry-run --headless --console")
- pipe_cmds.append(f"sed 's/{launcher_opts_find}/{launcher_opts_replace}'/")
- pipe_cmds.append(
- rf"sed 's/-Xmx[0-9a-z]\+/-Xmx{resources.mem_mb}m -Xms{resources.mem_mb}m/'"
- )
- pipe_cmds.append("tr --delete '\\n'")
- return "|".join(pipe_cmds) + f" > {output.launcher} && chmod a+x {output.launcher} "
-
-
-def get_macro_args_bigstitcher(wildcards, input, output):
- return "{dataset_xml} {pairwise_method} {ds_x} {ds_y} {ds_z} {do_filter} {min_r} {do_global} {global_strategy}".format(
- dataset_xml=output.dataset_xml,
- pairwise_method=config["bigstitcher"]["calc_pairwise_shifts"]["methods"][
- config["bigstitcher"]["calc_pairwise_shifts"]["method"]
- ],
- ds_x=config["bigstitcher"]["calc_pairwise_shifts"]["downsample_in_x"],
- ds_y=config["bigstitcher"]["calc_pairwise_shifts"]["downsample_in_y"],
- ds_z=config["bigstitcher"]["calc_pairwise_shifts"]["downsample_in_z"],
- do_filter=config["bigstitcher"]["filter_pairwise_shifts"]["enabled"],
- min_r=config["bigstitcher"]["filter_pairwise_shifts"]["min_r"],
- do_global=config["bigstitcher"]["global_optimization"]["enabled"],
- global_strategy=config["bigstitcher"]["global_optimization"]["strategies"][
- config["bigstitcher"]["global_optimization"]["strategy"]
- ],
- )
-
-
-def get_macro_args_zarr_fusion(wildcards, input, output):
- return "{dataset_xml} {downsampling} {channel:02d} {output_zarr} {bsx} {bsy} {bsz} {bsfx} {bsfy} {bsfz}".format(
- dataset_xml=input.dataset_xml,
- downsampling=config["bigstitcher"]["fuse_dataset"]["downsampling"],
- channel=get_stains(wildcards).index(wildcards.stain),
- output_zarr=output.zarr,
- bsx=config["bigstitcher"]["fuse_dataset"]["block_size_x"],
- bsy=config["bigstitcher"]["fuse_dataset"]["block_size_y"],
- bsz=config["bigstitcher"]["fuse_dataset"]["block_size_z"],
- bsfx=config["bigstitcher"]["fuse_dataset"]["block_size_factor_x"],
- bsfy=config["bigstitcher"]["fuse_dataset"]["block_size_factor_y"],
- bsfz=config["bigstitcher"]["fuse_dataset"]["block_size_factor_z"],
- )
-
def get_output_ome_zarr_uri():
if is_remote(config["root"]):
diff --git a/workflow/rules/import.smk b/workflow/rules/import.smk
index 8a9074a..421ff80 100644
--- a/workflow/rules/import.smk
+++ b/workflow/rules/import.smk
@@ -79,11 +79,6 @@ rule blaze_to_metadata_gcs:
rule blaze_to_metadata:
input:
ome_dir=get_input_dataset,
- params:
- in_tif_pattern=lambda wildcards, input: os.path.join(
- input.ome_dir,
- config["import_blaze"]["raw_tif_pattern"],
- ),
output:
metadata_json=temp(
bids(
@@ -197,10 +192,6 @@ rule tif_to_zarr:
ome_dir=get_input_dataset,
metadata_json=rules.copy_blaze_metadata.output.metadata_json,
params:
- in_tif_pattern=lambda wildcards, input: os.path.join(
- input.ome_dir,
- config["import_blaze"]["raw_tif_pattern"],
- ),
intensity_rescaling=config["import_blaze"]["intensity_rescaling"],
output:
zarr=temp(
@@ -289,7 +280,7 @@ rule tif_to_zarr_gcs:
),
group:
"preproc"
- threads: config["cores_per_rule"]
+ threads: 16
container:
config["containers"]["spimprep"]
script:
diff --git a/workflow/scripts/blaze_to_metadata.py b/workflow/scripts/blaze_to_metadata.py
index 608158e..25a27e3 100644
--- a/workflow/scripts/blaze_to_metadata.py
+++ b/workflow/scripts/blaze_to_metadata.py
@@ -4,25 +4,67 @@
import re
from itertools import product
from snakemake.io import glob_wildcards
-
-in_tif_pattern = snakemake.params.in_tif_pattern
+from glob import glob
+from pathlib import Path
+import os
+
+tif_files = glob(f"{snakemake.input.ome_dir}/*.tif")
+
+#check first tif file to see if it is zstack or not:
+if 'xyz-Table Z' in Path(tif_files[0]).name:
+ is_zstack=False
+else:
+ is_zstack=True
+
+#now check to see if there is just single tile
+if ' x ' in Path(tif_files[0]).name:
+ is_tiled=True
+else:
+ is_tiled=False
+
+if is_tiled:
+ if is_zstack:
+ in_tif_pattern = os.path.join(snakemake.input.ome_dir,snakemake.config["import_blaze"]["raw_tif_pattern_zstack"])
+ else:
+ in_tif_pattern = os.path.join(snakemake.input.ome_dir,snakemake.config["import_blaze"]["raw_tif_pattern"])
+else:
+ in_tif_pattern = os.path.join(snakemake.input.ome_dir,snakemake.config["import_blaze"]["raw_tif_pattern_notile"])
#add a wildcard constraint to ensure no
#subfolders get parsed (ie don't match anything with / in it):
prefix_constraint=r'[^/]+'
in_tif_pattern_constrained = in_tif_pattern.replace('{prefix}',f'{{prefix,{prefix_constraint}}}')
+
#parse the filenames to get number of channels, tiles etc..
-prefix, tilex, tiley, channel, zslice = glob_wildcards(in_tif_pattern_constrained)
+if is_tiled:
+ if is_zstack:
+ prefix, tilex, tiley, channel = glob_wildcards(in_tif_pattern_constrained,tif_files)
+ else:
+ prefix, tilex, tiley, channel, zslice = glob_wildcards(in_tif_pattern_constrained,tif_files)
+else:
+ prefix, channel, zslice = glob_wildcards(in_tif_pattern_constrained,tif_files)
+
+if is_tiled:
+ tiles_x = sorted(list(set(tilex)))
+ tiles_y = sorted(list(set(tiley)))
-tiles_x = sorted(list(set(tilex)))
-tiles_y = sorted(list(set(tiley)))
channels = sorted(list(set(channel)))
-zslices = sorted(list(set(zslice)))
prefixes = sorted(list(set(prefix)))
-#read in series metadata from first file
-in_tif = in_tif_pattern.format(tilex=tiles_x[0],tiley=tiles_y[0],prefix=prefixes[0],channel=channels[0],zslice=zslices[0])
+
+if not is_zstack:
+ zslices = sorted(list(set(zslice)))
+
+if is_tiled:
+ if is_zstack:
+ in_tif = in_tif_pattern.format(tilex=tiles_x[0],tiley=tiles_y[0],prefix=prefixes[0],channel=channels[0])
+ else:
+ #read in series metadata from first file
+ in_tif = in_tif_pattern.format(tilex=tiles_x[0],tiley=tiles_y[0],prefix=prefixes[0],channel=channels[0],zslice=zslices[0])
+else:
+ in_tif = in_tif_pattern.format(prefix=prefixes[0],channel=channels[0],zslice=zslices[0])
+
raw_tif = tifffile.TiffFile(in_tif,mode='r')
@@ -37,57 +79,81 @@
custom_metadata = ome_dict['OME']['Image']['ca:CustomAttributes']
-
#read tile configuration from the microscope metadata
if axes == 'CZYX':
- tile_config_pattern=r"Blaze\[(?P