register-pheno-pheno-qc

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,

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")

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,6 +199,23 @@ def test_ops_wdl(tmp_path):
     assert len(merge_sbs_metadata_df) > len(
         merge_sbs_metadata_df.query("~barcode_count_0.isna()")
     )
+    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",

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)

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])==str(reference_time):
+                        reference_time_index = i
+                        break
+            if reference_time_index == -1:
+                raise ValueError(f"Reference time {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 ''} \

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"
@@ -178,7 +183,8 @@ 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"
         String merge_meta_suffix = "merge-sbs-metadata"
@@ -206,7 +212,9 @@ 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
     String cell_intersects_boundary_directory_non_reference_t = output_stripped + cell_intersects_boundary_non_reference_t_suffix
 
@@ -427,7 +435,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=select_first([phenotype_dapi_channel_before_registration, 0]),
+                        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 +620,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,