diff --git a/wdl/starfish.wdl b/wdl/starfish.wdl deleted file mode 100644 index 23a5870ec..000000000 --- a/wdl/starfish.wdl +++ /dev/null @@ -1,178 +0,0 @@ -task register { - File input_tarball - Int upsampling - - command <<< - mkdir inputs - tar xf '${input_tarball}' -C inputs --strip-components 1 - - mkdir registered - - starfish register inputs/org.json registered --u '${upsampling}' - >>> - - runtime { - docker: "marcusczi/starfish:latest" - } - - output { - Array[File] registered_collection = glob("registered/*") - } -} - -task filter { - - Array[File] input_collection - Int disk_size - - command <<< - mkdir inputs - for input_file in ${sep=' ' input_collection}; do - mv $input_file inputs/ - done - - mkdir filtered - - starfish filter inputs/org.json filtered --ds '${disk_size}' - >>> - - runtime { - docker: "marcusczi/starfish:latest" - } - - output { - Array[File] filtered_collection = glob("filtered/*") - } -} - -task detect_spots { - - Array[File] input_collection - Int min_sigma - Int max_sigma - Int num_sigma - Float threshold - - command <<< - mkdir inputs - for input_file in ${sep=' ' input_collection}; do - mv $input_file inputs/ - done - - mkdir detected - - starfish detect_spots inputs/org.json detected dots --min_sigma '${min_sigma}' \ - --max_sigma '${max_sigma}' --num_sigma '${num_sigma}' --t '${threshold}' - >>> - - runtime { - docker: "marcusczi/starfish:latest" - } - - output { - File spots_geo = "detected/spots_geo.csv" - File encoder_table = "detected/encoder_table.csv" - } -} - -task segment { - - Array[File] input_collection - File spots_geo - File encoder_table - Float dapi_threshold - Float stain_threshold - Int minimum_distance - - command <<< - mkdir inputs - for input_file in ${sep=' ' input_collection}; do - mv $input_file inputs/ - done - - mkdir segmented - mv '${spots_geo}' segmented - mv '${encoder_table}' segmented - - - starfish segment inputs/org.json segmented stain --dt '${dapi_threshold}' \ - --st '${stain_threshold}' --md '${minimum_distance}' - >>> - - runtime { - docker: "marcusczi/starfish:latest" - } - - output { - Array[File] segmented_collection = glob("segmented/*") - } - -} - -task decode { - Array[File] input_collection - String decoder_type - - command <<< - mkdir inputs - for input_file in ${sep=' ' input_collection}; do - mv $input_file inputs/ - done - - starfish decode inputs --decoder_type '${decoder_type}' - >>> - - runtime { - docker: "marcusczi/starfish:latest" - } - - output { - Array[File] output_collection = glob("inputs/*") - } -} - - -workflow starfish { - - Int upsampling - Int disk_size - Int min_sigma - Int max_sigma - Int num_sigma - Float threshold - Float dapi_threshold - Float stain_threshold - Int minimum_distance - String decoder_type - - File input_tarball - - call register { - input: input_tarball=input_tarball, upsampling=upsampling - } - - call filter { - input: input_collection=register.registered_collection, disk_size=disk_size - } - - call detect_spots { - input: input_collection=filter.filtered_collection, min_sigma=min_sigma, - max_sigma=max_sigma, num_sigma=num_sigma, threshold=threshold - } - - call segment { - input: input_collection=filter.filtered_collection, dapi_threshold=dapi_threshold, - stain_threshold=stain_threshold, minimum_distance=minimum_distance, - spots_geo=detect_spots.spots_geo, encoder_table=detect_spots.encoder_table - } - - call decode { - input: input_collection=segment.segmented_collection, decoder_type=decoder_type - } - - output { - File spots_geo = detect_spots.spots_geo - File encoder_table = detect_spots.encoder_table - Array[File] decoded_results = decode.output_collection - } -} diff --git a/wdl/starfish_example_inputs.json b/wdl/starfish_example_inputs.json deleted file mode 100644 index b01513976..000000000 --- a/wdl/starfish_example_inputs.json +++ /dev/null @@ -1,13 +0,0 @@ -{ - "starfish.min_sigma": 4, - "starfish.threshold": 0.01, - "starfish.input_tarball": "/tmp/starfish/formatted.tar.gz", - "starfish.num_sigma": 20, - "starfish.disk_size": 10, - "starfish.upsampling": 1000, - "starfish.max_sigma": 6, - "starfish.dapi_threshold": 0.16, - "starfish.stain_threshold": 0.22, - "starfish.minimum_distance": 57, - "starfish.decoder_type": "iss" -} diff --git a/workflows/wdl/iss_published/inputs.json b/workflows/wdl/iss_published/inputs.json new file mode 100644 index 000000000..e4897a630 --- /dev/null +++ b/workflows/wdl/iss_published/inputs.json @@ -0,0 +1,5 @@ +{ + "Starfish.experiment": "https://d2nhj9g34unfro.cloudfront.net/browse/formatted/iss/20190506/experiment.json", + "Starfish.num_fovs": 15, + "Starfish.recipe_file": "https://raw.githubusercontent.com/spacetx/starfish/master/workflows/wdl/iss_published/recipe.py" +} \ No newline at end of file diff --git a/workflows/wdl/iss_published/recipe.py b/workflows/wdl/iss_published/recipe.py new file mode 100644 index 000000000..9883b1c79 --- /dev/null +++ b/workflows/wdl/iss_published/recipe.py @@ -0,0 +1,66 @@ +import starfish +from starfish import FieldOfView +from starfish.image import Filter +from starfish.image import ApplyTransform, LearnTransform +from starfish.spots import DetectSpots +from starfish.types import Axes + + +def process_fov(field_num: int, experiment_str: str): + """Process a single field of view of ISS data + Parameters + ---------- + field_num : int + the field of view to process + experiment_str : int + path of experiment json file + + Returns + ------- + DecodedSpots : + tabular object containing the locations of detected spots. + """ + + fov_str: str = f"fov_{int(field_num):03d}" + + # load experiment + experiment = starfish.Experiment.from_json(experiment_str) + + print(f"loading fov: {fov_str}") + fov = experiment[fov_str] + + # note the structure of the 5D tensor containing the raw imaging data + imgs = fov.get_image(FieldOfView.PRIMARY_IMAGES) + dots = fov.get_image("dots") + nuclei = fov.get_image("nuclei") + + masking_radius = 15 + print("Filter WhiteTophat") + filt = Filter.WhiteTophat(masking_radius, is_volume=False) + + filtered_imgs = filt.run(imgs, verbose=True, in_place=False) + filt.run(dots, verbose=True, in_place=True) + filt.run(nuclei, verbose=True, in_place=True) + + print("Learning Transform") + learn_translation = LearnTransform.Translation(reference_stack=dots, axes=Axes.ROUND, upsampling=1000) + transforms_list = learn_translation.run(imgs.max_proj(Axes.CH, Axes.ZPLANE)) + + print("Applying transform") + warp = ApplyTransform.Warp() + registered_imgs = warp.run(filtered_imgs, transforms_list=transforms_list, in_place=False, verbose=True) + + print("Detecting") + p = DetectSpots.BlobDetector( + min_sigma=1, + max_sigma=10, + num_sigma=30, + threshold=0.01, + measurement_type='mean', + ) + + intensities = p.run(registered_imgs, blobs_image=dots, blobs_axes=(Axes.ROUND, Axes.ZPLANE)) + + decoded = experiment.codebook.decode_per_round_max(intensities) + df = decoded.to_decoded_spots() + return df \ No newline at end of file diff --git a/workflows/wdl/iss_published/test_inputs.json b/workflows/wdl/iss_published/test_inputs.json new file mode 100644 index 000000000..6e0c8f9f8 --- /dev/null +++ b/workflows/wdl/iss_published/test_inputs.json @@ -0,0 +1,5 @@ +{ + "Starfish.experiment": "https://d2nhj9g34unfro.cloudfront.net/browse/formatted/iss/20190506/experiment.json", + "Starfish.num_fovs": 2, + "Starfish.recipe_file": "https://raw.githubusercontent.com/spacetx/starfish/master/workflows/wdl/iss_published/recipe.py" +} \ No newline at end of file diff --git a/workflows/wdl/iss_spaceTX/inputs.json b/workflows/wdl/iss_spaceTX/inputs.json new file mode 100644 index 000000000..ff28a53be --- /dev/null +++ b/workflows/wdl/iss_spaceTX/inputs.json @@ -0,0 +1,5 @@ +{ + "Starfish.experiment": "https://d2nhj9g34unfro.cloudfront.net/xiaoyan_qian/ISS_human_HCA_07_MultiFOV/main_files/experiment.json", + "Starfish.num_fovs": 539, + "Starfish.recipe_file": "https://raw.githubusercontent.com/spacetx/starfish/master/workflows/wdl/iss_spaceTX/recipe.py" +} \ No newline at end of file diff --git a/workflows/wdl/iss_spaceTX/recipe.py b/workflows/wdl/iss_spaceTX/recipe.py new file mode 100644 index 000000000..f2f2630dd --- /dev/null +++ b/workflows/wdl/iss_spaceTX/recipe.py @@ -0,0 +1,55 @@ +import numpy as np +import starfish +from starfish.types import Axes + + +def process_fov(field_num: int, experiment_str: str): + """Process a single field of view of ISS data + Parameters + ---------- + field_num : int + the field of view to process + experiment_str : int + path of experiment json file + + Returns + ------- + DecodedSpots : + tabular object containing the locations of detected spots. + """ + fov_str: str = f"fov_{int(field_num):03d}" + + # load experiment + experiment = starfish.Experiment.from_json(experiment_str) + + fov = experiment[fov_str] + imgs = fov.get_image(starfish.FieldOfView.PRIMARY_IMAGES) + dots = imgs.max_proj(Axes.CH) + + # filter + filt = starfish.image.Filter.WhiteTophat(masking_radius=15, is_volume=False) + filtered_imgs = filt.run(imgs, verbose=True, in_place=False) + filt.run(dots, verbose=True, in_place=True) + + # find threshold + tmp = dots.sel({Axes.ROUND:0, Axes.CH:0, Axes.ZPLANE:0}) + dots_threshold = np.percentile(np.ravel(tmp.xarray.values), 50) + + # find spots + p = starfish.spots.DetectSpots.BlobDetector( + min_sigma=1, + max_sigma=10, + num_sigma=30, + threshold=dots_threshold, + measurement_type='mean', + ) + + # blobs = dots; define the spots in the dots image, but then find them again in the stack. + intensities = p.run(filtered_imgs, blobs_image=dots, blobs_axes=(Axes.ROUND, Axes.ZPLANE)) + + # decode + decoded = experiment.codebook.decode_per_round_max(intensities) + + # save results + df = decoded.to_decoded_spots() + return df \ No newline at end of file diff --git a/workflows/wdl/iss_spaceTX/test_inputs.json b/workflows/wdl/iss_spaceTX/test_inputs.json new file mode 100644 index 000000000..2ad48b553 --- /dev/null +++ b/workflows/wdl/iss_spaceTX/test_inputs.json @@ -0,0 +1,5 @@ +{ + "Starfish.experiment": "https://d2nhj9g34unfro.cloudfront.net/xiaoyan_qian/ISS_human_HCA_07_MultiFOV/main_files/experiment.json", + "Starfish.num_fovs": 2, + "Starfish.recipe_file": "https://raw.githubusercontent.com/spacetx/starfish/master/workflows/wdl/iss_spaceTX/recipe.py" +} \ No newline at end of file diff --git a/workflows/wdl/merfish_published/inputs.json b/workflows/wdl/merfish_published/inputs.json new file mode 100644 index 000000000..527427274 --- /dev/null +++ b/workflows/wdl/merfish_published/inputs.json @@ -0,0 +1,5 @@ +{ + "Starfish.experiment": "https://d2nhj9g34unfro.cloudfront.net/browse/formatted/MERFISH/20190511/experiment.json", + "Starfish.num_fovs": 495, + "Starfish.recipe_file": "https://raw.githubusercontent.com/spacetx/starfish/saxelrod-standard-wdl/workflows/wdl/merfish_published/recipe.py" +} \ No newline at end of file diff --git a/workflows/wdl/merfish_published/recipe.py b/workflows/wdl/merfish_published/recipe.py new file mode 100644 index 000000000..b027ebb10 --- /dev/null +++ b/workflows/wdl/merfish_published/recipe.py @@ -0,0 +1,72 @@ +import numpy as np +from copy import deepcopy + +import starfish +from starfish import FieldOfView +from starfish.types import Features, Axes +from starfish.image import Filter +from starfish.types import Clip +from starfish.spots import DetectPixels + + +def process_fov(field_num: int, experiment_str: str): + """Process a single field of view of MERFISH data + Parameters + ---------- + field_num : int + the field of view to process + experiment_str : int + path of experiment json file + + Returns + ------- + DecodedSpots : + tabular object containing the locations of detected spots. + """ + fov_str: str = f"fov_{int(field_num):03d}" + # load experiment + experiment = starfish.Experiment.from_json(experiment_str) + + print(f"Loading fov: {fov_str}") + fov = experiment[fov_str] + imgs = fov.get_image(FieldOfView.PRIMARY_IMAGES) + + print("Gaussian High Pass") + ghp = Filter.GaussianHighPass(sigma=3) + high_passed = ghp.run(imgs, verbose=True, in_place=False) + + print("Deconvolve") + dpsf = Filter.DeconvolvePSF(num_iter=15, sigma=2, clip_method=Clip.SCALE_BY_CHUNK) + deconvolved = dpsf.run(high_passed, verbose=True, in_place=False) + + print("Guassian Low Pass") + glp = Filter.GaussianLowPass(sigma=1) + low_passed = glp.run(deconvolved, in_place=False, verbose=True) + + scale_factors = { + (t[Axes.ROUND], t[Axes.CH]): t['scale_factor'] + for t in experiment.extras['scale_factors'] + } + filtered_imgs = deepcopy(low_passed) + + for selector in imgs._iter_axes(): + data = filtered_imgs.get_slice(selector)[0] + scaled = data / scale_factors[selector[Axes.ROUND.value], selector[Axes.CH.value]] + filtered_imgs.set_slice(selector, scaled, [Axes.ZPLANE]) + + print("Decode") + psd = DetectPixels.PixelSpotDecoder( + codebook=experiment.codebook, + metric='euclidean', # distance metric to use for computing distance between a pixel vector and a codeword + norm_order=2, # the L_n norm is taken of each pixel vector and codeword before computing the distance. this is n + distance_threshold=0.5176, # minimum distance between a pixel vector and a codeword for it to be called as a gene + magnitude_threshold=1.77e-5, # discard any pixel vectors below this magnitude + min_area=2, # do not call a 'spot' if it's area is below this threshold (measured in pixels) + max_area=np.inf, # do not call a 'spot' if it's area is above this threshold (measured in pixels) + ) + + initial_spot_intensities, prop_results = psd.run(filtered_imgs) + + spot_intensities = initial_spot_intensities.loc[initial_spot_intensities[Features.PASSES_THRESHOLDS]] + df = spot_intensities.to_decoded_spots() + return df \ No newline at end of file diff --git a/workflows/wdl/merfish_published/test_inputs.json b/workflows/wdl/merfish_published/test_inputs.json new file mode 100644 index 000000000..1bf188ee4 --- /dev/null +++ b/workflows/wdl/merfish_published/test_inputs.json @@ -0,0 +1,5 @@ +{ + "Starfish.experiment": "https://d2nhj9g34unfro.cloudfront.net/browse/formatted/MERFISH/20190511/experiment.json", + "Starfish.num_fovs": 2, + "Starfish.recipe_file": "https://raw.githubusercontent.com/spacetx/starfish/saxelrod-standard-wdl/workflows/wdl/merfish_published/recipe.py" +} \ No newline at end of file diff --git a/workflows/wdl/starfish.wdl b/workflows/wdl/starfish.wdl new file mode 100644 index 000000000..0dcd8a91c --- /dev/null +++ b/workflows/wdl/starfish.wdl @@ -0,0 +1,100 @@ + +task process_field_of_view { + + String experiment + Int field_of_view + String recipe_file + + command <<< + wget -O recipe.py ${recipe_file} + + python3 <>> + + runtime { + docker: "spacetx/starfish:0.1.0-simple" + memory: "16 GB" + cpu: "4" + disk: "local-disk 100 SDD" + } + + output { + File decoded_csv = "decoded.csv" + } +} + + +task concatenate_fovs { + Array[File] decoded_csvs + + command <<< + python <>> + + runtime { + docker: "spacetx/starfish:0.1.0-simple" + memory: "16 GB" + cpu: "1" + disk: "local-disk 100 SDD" + } + + output { + File decoded_spots = "decoded_concatenated.csv" + } +} + + +workflow Starfish { + Int num_fovs + String experiment + String recipe_file + + Array[Int] fields_of_view = range(num_fovs) + + scatter(fov in fields_of_view) { + call process_field_of_view { + input: + experiment = experiment, + field_of_view = fov, + recipe_file = recipe_file + } + } + + call concatenate_fovs { + input: + decoded_csvs = process_field_of_view.decoded_csv + } + + output { + File decoded_spots = concatenate_fovs.decoded_spots + } +}