From 0c9e242fb0193049d33413e169044b43cd8d2f0b Mon Sep 17 00:00:00 2001 From: Joshua Gould Date: Mon, 13 Jul 2026 13:52:12 -0400 Subject: [PATCH 1/4] register pheno to pheno qc --- scallops/registration/itk.py | 22 +++++++++----- scallops/tests/test_wdl.py | 15 ++++++---- scallops/utils.py | 55 +++++++++++++--------------------- wdl/ops_tasks.wdl | 57 ++++++++++++++++++++++++++++-------- wdl/ops_workflow.wdl | 31 ++++++++++++++++++++ 5 files changed, 119 insertions(+), 61 deletions(-) diff --git a/scallops/registration/itk.py b/scallops/registration/itk.py index b36501d..7bad5ab 100644 --- a/scallops/registration/itk.py +++ b/scallops/registration/itk.py @@ -596,7 +596,7 @@ def _write_value(x, idx, val): dim="c", keep_attrs=True ) if isinstance(moving_image, xr.DataArray): - times = moving_image.coords["t"].values + times = moving_image.coords["t"].values.tolist() else: times = [] for img in moving_image: @@ -621,11 +621,17 @@ def _write_value(x, idx, val): ) channels = ( - [moving_image.coords["c"].values for i in range(moving_image.sizes["t"])] + [ + moving_image.coords["c"].values.tolist() + for _ in range(moving_image.sizes["t"]) + ] if isinstance(moving_image, xr.DataArray) - else [moving_image[i].coords["c"].values for i in range(len(moving_image))] + else [ + moving_image[i].coords["c"].values.tolist() + for i in range(len(moving_image)) + ] ) - + result_image_attrs = moving_image_reference_t_moving_c.attrs.copy() if unroll_channels: logger.info( f"Channels per time point: {', '.join([str(c) for c in nchannels])}." @@ -646,6 +652,9 @@ def _write_value(x, idx, val): all_channels.append(channel_name) coords = dict(c=all_channels) leading_dims_shape = (len(all_channels),) + result_image_attrs["src_channels"] = channels + result_image_attrs["src_times"] = times + else: # all images must have the same number of channels dims = ["t", "c", "y", "x"] leading_dims_shape = len(times), nchannels[0] @@ -662,7 +671,7 @@ def _write_value(x, idx, val): coords=coords, chunk_size=_get_chunk_size(moving_image_reference_t_moving_c), shape=result_image_shape, - attrs=moving_image_reference_t_moving_c.attrs.copy(), + attrs=result_image_attrs, dtype=moving_image_reference_t_moving_c.dtype, ) ) @@ -1151,7 +1160,6 @@ def _itk_transform_image_zarr( chunksize: tuple[int, int] | None, image_attrs: dict, image_spacing: None | tuple[float, float] = None, - channels_transform_parameter_objects: dict[int, itk.ParameterObject] | None = None, ): """Apply an ITK-based transformation to an image and store the result in a Zarr group. @@ -1162,8 +1170,6 @@ def _itk_transform_image_zarr( :param chunksize: Chunk size to use for the Zarr dataset (default is (1024, 1024)). :param image_attrs: Attributes to store with the transformed image. :param image_spacing: Physical spacing of the input image. - :param channels_transform_parameter_objects: Maps channel index to ITK parameter - object """ if not isinstance(image_root, zarr.Group): image_root = open_ome_zarr(image_root, mode="a") diff --git a/scallops/tests/test_wdl.py b/scallops/tests/test_wdl.py index 0e1b2ee..ad54674 100644 --- a/scallops/tests/test_wdl.py +++ b/scallops/tests/test_wdl.py @@ -141,12 +141,12 @@ def test_ops_wdl(tmp_path): ) phenotype_tile[10, 10] = 2 exp = Experiment( - images={"A1-102-1": pheno_img, "A1-102-2": pheno_img}, + images={"A1-102-IF": pheno_img, "A1-102-FISH": pheno_img}, labels={ - "A1-102-1-mask": phenotype_mask, - "A1-102-1-tile": phenotype_tile, - "A1-102-2-mask": phenotype_mask, - "A1-102-2-tile": phenotype_tile, + "A1-102-IF-mask": phenotype_mask, + "A1-102-IF-tile": phenotype_tile, + "A1-102-FISH-mask": phenotype_mask, + "A1-102-FISH-tile": phenotype_tile, }, ) exp.save(str(pheno_dir)) @@ -161,6 +161,7 @@ def test_ops_wdl(tmp_path): "pheno_registration_extra_arguments": "--no-landmarks", "phenotype_cyto_channel": [1], "phenotype_dapi_channel": 0, + "reference_phenotype_time": "IF", "phenotype_url": str(pheno_dir.absolute()), "phenotype_nuclei_features": ["intensity_0", "intensity_1"], # 2 batches @@ -198,7 +199,9 @@ def test_ops_wdl(tmp_path): assert len(merge_sbs_metadata_df) > len( merge_sbs_metadata_df.query("~barcode_count_0.isna()") ) - + # assert merge_sbs_metadata_df.columns.str.contains('qc').sum() == xxx + # assert merge_sbs_metadata_df.columns.str.contains('qc').sum() == xxx + # assert merge_sbs_metadata_df.columns.str.contains('qc').sum() == xxx for col in [ "Nuclei_AreaShape_Area", "Cells_AreaShape_Area", diff --git a/scallops/utils.py b/scallops/utils.py index f015589..9dc1f85 100644 --- a/scallops/utils.py +++ b/scallops/utils.py @@ -688,7 +688,7 @@ def _list_images_wdl( """Used by WDL workflow to output info about images""" from scallops.io import _set_up_experiment - batch_size = 0 + batch_size = 1 expected_cycles = None if expected_cycles_str != "": expected_cycles = int(expected_cycles_str) @@ -710,41 +710,26 @@ def _list_images_wdl( if not save_group_size: with open("group_size.txt", "wt") as f: f.write("0\n") - if batch_size > 0: - with open("groups.txt", "wt") as f: - ids = [] - first = True - for g, file_list, metadata in exp_gen: - if first: - first = False - if save_group_size: - _write_group_size(metadata) - if not groupby_t and "t" in metadata["file_metadata"][0]: - t = [md["t"] for md in metadata["file_metadata"]] - if expected_cycles is not None: - assert len(t) == expected_cycles - - ids.append('"' + metadata["id"] + '"') - if len(ids) == batch_size: - f.write(" ".join(ids)) - f.write("\n") - ids = [] - if len(ids) > 0: - f.write(" ".join(ids)) - f.write("\n") - else: - with open("groups.txt", "wt") as f: - first = True - for g, file_list, metadata in exp_gen: - f.write(metadata["id"]) - f.write("\n") - if first: - first = False - if save_group_size: - _write_group_size(metadata) - if not groupby_t and "t" in metadata["file_metadata"][0]: - t = [md["t"] for md in metadata["file_metadata"]] + first = True + subset_ids = [] + for g, file_list, metadata in exp_gen: + subset_ids.append(metadata["id"]) + + if first: + first = False + if save_group_size: + _write_group_size(metadata) + if not groupby_t and "t" in metadata["file_metadata"][0]: + t = [md["t"] for md in metadata["file_metadata"]] + if expected_cycles is not None: + assert len(t) == expected_cycles + + with open("groups.txt", "wt") as f: # ["plate1-A1", "plate1-A2", ...] + for i in range(0, len(subset_ids), batch_size): + selected = subset_ids[i : i + batch_size] + f.write(" ".join(selected)) + f.write("\n") with open("groupby.txt", "wt") as f: for g in groupby: f.write(g) diff --git a/wdl/ops_tasks.wdl b/wdl/ops_tasks.wdl index f17dc52..6eb9b64 100644 --- a/wdl/ops_tasks.wdl +++ b/wdl/ops_tasks.wdl @@ -275,9 +275,10 @@ task register_qc { Int channel String subset String output_directory - String channel_prefix + String? channel_prefix Array[String] groupby Boolean? force + String? reference_time String docker String zones @@ -306,8 +307,8 @@ task register_qc { image_pattern = "~{image_pattern}" label_type = "~{label_type}" labels = "~{labels}" - - channel = "~{channel}" + reference_time = "~{reference_time}" + channel = ~{channel} subset = "~{subset}".split(" ") output_directory = "~{output_directory}" groupby = "~{sep=',' groupby}".split(",") @@ -319,16 +320,45 @@ task register_qc { if len(keys) == 0: raise ValueError("No images found") image = exp.images[keys[0]] - size_t = image.sizes['t'] - size_c = image.sizes['c'] - channel_rename = {} - t = 0 - for i in range(int(channel), size_t * size_c, size_c): - channel_rename[f"{i}"] = f"{channel_prefix}{t}" - t += 1 + channel_rename = None + if "src_channels" in image.attrs: + src_channels = image.attrs["src_channels"] + src_times = image.attrs["src_times"] + ref_channel = channel + reference_time_index = 0 + if reference_time != "": + reference_time_index = -1 + for i in range(len(src_times)): + if str(src_times[i])==reference_time: + reference_time_index = -1 + break + if reference_time_index == -1: + raise ValueError("Reference time not found") + test_channels = [] + channel_index = channel + channel_rename = {} + + for i in range(len(src_times)): + n_channels = len(src_channels[i]) + if reference_time_index == i: + ref_channel = str(channel_index) + else: + test_channels.append(str(channel_index)) + channel_rename[f"{channel_index}"] = f"{src_times[i]}" + channel_index += n_channels + features = f"correlationpearsonbox_{ref_channel}_{','.join(test_channels)}" + else: + size_t = image.sizes['t'] + size_c = image.sizes['c'] + channel_rename = {} + t = 0 + for i in range(channel, size_t * size_c, size_c): + channel_rename[f"{i}"] = f"{channel_prefix}{t}" + t += 1 + features = f"correlationpearsonbox_{channel}_{channel}:{size_t * size_c}:{size_c}" cmd = ["scallops", "features"] - cmd += [f"--features-{label_type}", f"correlationpearsonbox_{channel}_{channel}:{size_t * size_c}:{size_c}"] + cmd += [f"--features-{label_type}", features] cmd += ["--labels", labels] if image_pattern != "": @@ -340,7 +370,8 @@ task register_qc { cmd += subset cmd += ["--output", output_directory] cmd += ["--images", images] - cmd += ["--channel-rename", f"{json.dumps(channel_rename)}"] + if channel_rename is not None: + cmd += ["--channel-rename", f"{json.dumps(channel_rename)}"] if force == "true": cmd.append("--force") @@ -697,6 +728,7 @@ task merge { String? objects_cell String? register_pheno_to_iss_qc String? register_iss_to_iss_qc + String? register_pheno_to_pheno_qc String? objects_cytosol String? cell_intersects_boundary String? cell_intersects_boundary_t @@ -737,6 +769,7 @@ task merge { ~{cell_intersects_boundary_t} \ ~{register_pheno_to_iss_qc} \ ~{register_iss_to_iss_qc} \ + ~{register_pheno_to_pheno_qc} \ --subset ~{subset} \ ~{"--barcode-col " + barcode_column} \ ~{if defined(extra_arguments) then extra_arguments else ''} \ diff --git a/wdl/ops_workflow.wdl b/wdl/ops_workflow.wdl index 1955d4a..e57b0bf 100644 --- a/wdl/ops_workflow.wdl +++ b/wdl/ops_workflow.wdl @@ -98,6 +98,7 @@ workflow ops_workflow { Boolean force_find_objects = false Boolean force_register_pheno_to_iss_qc = false Boolean force_register_iss_to_iss_qc = false + Boolean force_register_pheno_to_pheno_qc = false # general options Array[String]? subset @@ -144,6 +145,10 @@ workflow ops_workflow { Int register_pheno_to_iss_qc_cpu = 48 String register_pheno_to_iss_qc_disks = "local-disk 200 HDD" + String register_pheno_to_pheno_qc_memory = "96 GiB" + Int register_pheno_to_pheno_qc_cpu = 48 + String register_pheno_to_pheno_qc_disks = "local-disk 200 HDD" + String register_iss_to_iss_qc_memory = "48 GiB" Int register_iss_to_iss_qc_cpu = 24 String register_iss_to_iss_qc_disks = "local-disk 200 HDD" @@ -179,6 +184,7 @@ workflow ops_workflow { String register_pheno_to_pheno_transform_suffix = "pheno-to-pheno-transforms" String register_pheno_to_iss_qc_suffix = "pheno-to-iss-qc" String register_iss_to_iss_qc_directory = "iss-to-iss-qc" + String register_pheno_to_pheno_qc_suffix = "pheno-to-pheno-qc" String spot_detect_suffix = "spot-detect.zarr" String reads_suffix = "reads" String merge_meta_suffix = "merge-sbs-metadata" @@ -207,6 +213,7 @@ workflow ops_workflow { String merge_meta_directory = output_stripped + merge_meta_suffix String merge_features_directory = output_stripped + merge_features_suffix String register_pheno_to_iss_qc_directory = output_stripped + register_pheno_to_iss_qc_suffix + String register_pheno_to_pheno_qc_directory = output_stripped + register_pheno_to_pheno_qc_suffix String cell_intersects_boundary_directory = output_stripped + cell_intersects_boundary_suffix String cell_intersects_boundary_directory_non_reference_t = output_stripped + cell_intersects_boundary_non_reference_t_suffix @@ -427,7 +434,30 @@ workflow ops_workflow { max_retries = max_retries } } + if(length(times)>1) { + call tasks.register_qc as register_pheno_to_pheno_qc { + input: + images=select_first([register_pheno_to_pheno.moving_output_url]), + image_pattern=image_pattern_after_registration, + channel=phenotype_dapi_channel, + label_type='nuclei', + reference_time=reference_phenotype_time, + output_directory=register_pheno_to_pheno_qc_directory, + labels=select_first([segment_nuclei.output_url]), + subset = group, + groupby=groupby, + force = force_register_pheno_to_pheno_qc, + docker=docker, + zones = zones, + preemptible = preemptible, + aws_queue_arn = aws_queue_arn, + disks = register_pheno_to_pheno_qc_disks, + memory = register_pheno_to_pheno_qc_memory, + cpu = register_pheno_to_pheno_qc_cpu, + max_retries = max_retries + } + } if(iss_url_supplied && pheno_url_supplied) { call tasks.register_elastix as register_pheno_to_iss { input: @@ -589,6 +619,7 @@ workflow ops_workflow { cell_intersects_boundary=cell_intersects_boundary.output_url, cell_intersects_boundary_t=cell_intersects_boundary_t.output_url, register_pheno_to_iss_qc=register_pheno_to_iss_qc.output_url, + register_pheno_to_pheno_qc=register_pheno_to_pheno_qc.output_url, register_iss_to_iss_qc=register_iss_to_iss_qc.output_url, barcodes=select_first([barcodes]), barcode_column=barcode_column, From bedb7cb3aad544ee134a0f30c1ce52303ba78db7 Mon Sep 17 00:00:00 2001 From: Joshua Gould Date: Mon, 13 Jul 2026 14:05:44 -0400 Subject: [PATCH 2/4] register pheno to pheno qc --- wdl/ops_tasks.wdl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/wdl/ops_tasks.wdl b/wdl/ops_tasks.wdl index 6eb9b64..73cb69a 100644 --- a/wdl/ops_tasks.wdl +++ b/wdl/ops_tasks.wdl @@ -329,11 +329,11 @@ task register_qc { if reference_time != "": reference_time_index = -1 for i in range(len(src_times)): - if str(src_times[i])==reference_time: - reference_time_index = -1 + if str(src_times[i])==str(reference_time): + reference_time_index = i break if reference_time_index == -1: - raise ValueError("Reference time not found") + raise ValueError(f"Reference time {reference_time} not found.") test_channels = [] channel_index = channel channel_rename = {} From b518e6639200598e2822a8cbbd1142690911634b Mon Sep 17 00:00:00 2001 From: Joshua Gould Date: Mon, 13 Jul 2026 14:26:56 -0400 Subject: [PATCH 3/4] register pheno to pheno qc --- scallops/tests/test_wdl.py | 21 ++++++++++++++++++--- wdl/ops_workflow.wdl | 5 +++-- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/scallops/tests/test_wdl.py b/scallops/tests/test_wdl.py index ad54674..57dbc3b 100644 --- a/scallops/tests/test_wdl.py +++ b/scallops/tests/test_wdl.py @@ -199,9 +199,24 @@ def test_ops_wdl(tmp_path): assert len(merge_sbs_metadata_df) > len( merge_sbs_metadata_df.query("~barcode_count_0.isna()") ) - # assert merge_sbs_metadata_df.columns.str.contains('qc').sum() == xxx - # assert merge_sbs_metadata_df.columns.str.contains('qc').sum() == xxx - # assert merge_sbs_metadata_df.columns.str.contains('qc').sum() == xxx + bbox_cols = sorted( + merge_sbs_metadata_df.columns[ + merge_sbs_metadata_df.columns.str.contains("PearsonBox") + ].tolist() + ) + assert bbox_cols == [ + "Nuclei_Correlation_PearsonBox_FISH_IF", + "Nuclei_Correlation_PearsonBox_ISS0_ISS1", + "Nuclei_Correlation_PearsonBox_ISS0_ISS2", + "Nuclei_Correlation_PearsonBox_ISS0_ISS3", + "Nuclei_Correlation_PearsonBox_ISS0_ISS4", + "Nuclei_Correlation_PearsonBox_ISS0_ISS5", + "Nuclei_Correlation_PearsonBox_ISS0_ISS6", + "Nuclei_Correlation_PearsonBox_ISS0_ISS7", + "Nuclei_Correlation_PearsonBox_ISS0_ISS8", + "Nuclei_Correlation_PearsonBox_ISS_PHENO", + ] + for col in [ "Nuclei_AreaShape_Area", "Cells_AreaShape_Area", diff --git a/wdl/ops_workflow.wdl b/wdl/ops_workflow.wdl index e57b0bf..c2f43c8 100644 --- a/wdl/ops_workflow.wdl +++ b/wdl/ops_workflow.wdl @@ -183,7 +183,7 @@ workflow ops_workflow { String register_pheno_to_pheno_suffix = "pheno-registered.zarr" String register_pheno_to_pheno_transform_suffix = "pheno-to-pheno-transforms" String register_pheno_to_iss_qc_suffix = "pheno-to-iss-qc" - String register_iss_to_iss_qc_directory = "iss-to-iss-qc" + String register_iss_to_iss_qc_suffix = "iss-to-iss-qc" String register_pheno_to_pheno_qc_suffix = "pheno-to-pheno-qc" String spot_detect_suffix = "spot-detect.zarr" String reads_suffix = "reads" @@ -212,6 +212,7 @@ workflow ops_workflow { String reads_directory = output_stripped + reads_suffix String merge_meta_directory = output_stripped + merge_meta_suffix String merge_features_directory = output_stripped + merge_features_suffix + String register_iss_to_iss_qc_directory = output_stripped + register_iss_to_iss_qc_suffix String register_pheno_to_iss_qc_directory = output_stripped + register_pheno_to_iss_qc_suffix String register_pheno_to_pheno_qc_directory = output_stripped + register_pheno_to_pheno_qc_suffix String cell_intersects_boundary_directory = output_stripped + cell_intersects_boundary_suffix @@ -440,7 +441,7 @@ workflow ops_workflow { input: images=select_first([register_pheno_to_pheno.moving_output_url]), image_pattern=image_pattern_after_registration, - channel=phenotype_dapi_channel, + channel=select_first([phenotype_dapi_channel_before_registration, 0]), label_type='nuclei', reference_time=reference_phenotype_time, output_directory=register_pheno_to_pheno_qc_directory, From 6563e97470a79417d0baa8de239e73ef17a65122 Mon Sep 17 00:00:00 2001 From: Joshua Gould Date: Mon, 13 Jul 2026 15:17:49 -0400 Subject: [PATCH 4/4] Fixed signature --- scallops/cli/register.py | 1 - 1 file changed, 1 deletion(-) diff --git a/scallops/cli/register.py b/scallops/cli/register.py index 9320482..706e7ad 100644 --- a/scallops/cli/register.py +++ b/scallops/cli/register.py @@ -526,7 +526,6 @@ def transform_all_images( _itk_transform_image_zarr( image=array, transform_parameter_object=transform_parameter_object, - channels_transform_parameter_objects=channels_transform_parameter_objects, image_attrs=dict(), image_spacing=moving_image_spacing, image_name=image_key,