ISS decoding module

config/defaults.yml

@@ -9,12 +9,14 @@ workflow:
   segmentation: unmicst
   segmentation-recyze: false
   downstream: scimap
+  iss_decoding: false
 options:
   ashlar: -m 30
   cypository: --model zeisscyto
   ilastik: --num_channels 1
   mcquant: --masks cell*.tif
   naivestates: -p png
+  starfish: --tile-size 2000
 modules:
   illumination:
     name: basic
@@ -114,3 +116,7 @@ modules:
     name: autominerva
     container: labsyspharm/auto-minerva
     version: '2022-06-06'
+  iss_decoding:
+    name: starfish
+    container: spacetx/starfish
+    version: 0.2.2

config/schema.yml

@@ -13,6 +13,7 @@ workflow:
   - downstream
   - ilastik-model
   - mesmer-model
+  - iss_decoding
 deprecated:
   quantification-mask: "--quant-opts '--masks ...'"
   illum: --start-at illumination

docs/parameters/other.md

@@ -14,10 +14,14 @@ Segmentation
 1. [Cypository](./other.html#cypository)
 1. [Mesmer](./other.html#mesmer)
 
-Clsutering and cell type inference
+Clustering and cell type inference
 1. [Clustering](./other.html#clustering) 
 1. [naivestates](./other.html#naivestates)
 
+In-situ sequencing spot deconvolution
+1. [Starfish](./other.html#starfish)
+
+
 <br>
 
 [Back to main modules](./){: .btn .btn-outline} 
@@ -258,3 +262,30 @@ Nextflow will write all outputs to the `cell-states/naivestates/` subdirectory w
 |`--mct <filename>` | |The tool has a basic marker -> cell type (mct) mapping in `typemap.csv`. More sophisticated mct mappings can be defined by creating a `custom-map.csv` file with two columns: `Marker` and `State`. |
 
 [Back to top](./other.html#other-modules){: .btn .btn-purple} [Back to main modules](./){: .btn .btn-outline} 
+
+## Starfish
+
+### Description
+starfish is a Python library for processing images of image-based spatial transcriptomics. We implemented part of the starfish pipeline to process ISS images starting from ... until the final result, a table with locations of all called spots.
+### Usage
+Add a `iss_decoding:` field to [workflow parameters]({{site.baseurl}}/parameters/) to select starfish.
+
+* Example `params.yml`:
+
+``` yaml
+workflow:
+  stop-at: iss_decoding
+  downstream: naivestates
+  naivestates-model: /full/path/to/mct.csv
+options:
+  naivestates: --log no
+```
+* Default naivestates options: `-p png`
+* Running outside of MCMICRO: [Instructions](https://github.com/labsyspharm/mcmicro-ilastik){:target="_blank"}.
+### Inputs
+
+
+### Outputs
+
+
+### Optional arguments

lib/mcmicro/Flow.groovy

@@ -32,7 +32,8 @@ static def flowSegment(wfp) {
         "segmentation",     // Step 4
         "watershed",        // Step 5
         "quantification",   // Step 6
-        "downstream"]       // Step 7
+        "downstream",       // Step 7
+        "iss_decoding"]   // Step 8
 
     // Identify starting and stopping indices
     int idxStart = mcsteps.indexOf( wfp['start-at'] )
@@ -65,7 +66,8 @@ static def precomputed(wfp) {
         dearray:            idxStart > 3 && wfp.tma,
         'probability-maps': idxStart == 5,
         segmentation:       idxStart == 6,
-        quantification:     idxStart == 7
+        quantification:     idxStart == 7,
+        iss_decoding:       idxStart == 8
     ]
 }
 
@@ -94,6 +96,8 @@ static def doirun(step, wfp) {
             return(idxStart <= 6 && idxStop >= 6)
         case 'downstream':
             return(idxStart <= 7 && idxStop >= 7)
+        case 'iss_decoding':
+            return(idxStart <= 8 && idxStop >= 8)
         case 'viz':
             return(wfp.viz)
         default:

main.nf

@@ -49,6 +49,9 @@ findFiles = { key, pattern, ife -> pre[key] ?
     Channel.fromPath("${params.in}/$key/$pattern").ifEmpty(ife) : Channel.empty()
 }
 
+// Find spot images in starfish_input
+starimgs = Channel.fromPath( "${params.in}/primary", checkIfExists: true )
+
 // Some image formats store multiple fields of view in a single file. Other
 // formats store each field separately, typically in .tif files, with a separate
 // index file to tie them together. We will look for the index files from
@@ -98,6 +101,7 @@ include {illumination}   from "$projectDir/modules/illumination"
 include {registration}   from "$projectDir/modules/registration"
 include {dearray}        from "$projectDir/modules/dearray"
 include {segmentation}   from "$projectDir/modules/segmentation"
+include {starfish}       from "$projectDir/modules/iss_decoding"
 include {quantification} from "$projectDir/modules/quantification"
 include {downstream}     from "$projectDir/modules/downstream"
 include {viz}            from "$projectDir/modules/viz"
@@ -124,6 +128,9 @@ workflow {
     tmacores = dearray.out.cores.mix(pre_cores)
     tmamasks = dearray.out.masks.mix(pre_masks)
 
+    // Is the data type ISS?
+    starfish(mcp, starimgs)
+
     // Reconcile WSI and TMA processing for downstream segmentation
     allimg = img.wsi.mix(tmacores)
     segmentation(mcp, allimg, tmamasks, pre_pmap)

modules/iss_decoding.nf

@@ -0,0 +1,138 @@
+// Import utility functions from lib/mcmicro/*.groovy
+import mcmicro.*
+
+process starfish_tile {
+
+    container "${params.contPfx}${module.container}:${module.version}"
+    publishDir "${params.in}/iss_processing", mode: 'copy', pattern: "*.tif"
+
+    publishDir "${Flow.QC(params.in, 'provenance')}", mode: 'copy', 
+      pattern: '.command.{sh,log}',
+      saveAs: {fn -> fn.replace('.command', "${module.name}-${task.index}")}
+
+  input:
+    val mcp
+    val module
+    path code_tile
+    path img_dir
+
+  output:
+    path 'TILED'
+    tuple path('.command.sh'), path('.command.log')
+
+  when: mcp.workflow["iss_decoding"]
+
+    """
+    python $code_tile --input ${img_dir} --output TILED
+    """
+}
+
+process starfish_convert {
+
+    container "${params.contPfx}${module.container}:${module.version}"
+    publishDir "${params.in}/iss_processing", mode: 'copy', pattern: "*.tif"
+
+    publishDir "${Flow.QC(params.in, 'provenance')}", mode: 'copy',
+      pattern: '.command.{sh,log}',
+      saveAs: {fn -> fn.replace('.command', "${module.name}-${task.index}")}
+
+  input:
+    val mcp
+    val module
+    path code
+    path TILED
+
+  output:
+    path 'SpaceTx'
+    tuple path('.command.sh'), path('.command.log')
+
+  when: mcp.workflow["iss_decoding"]
+
+    """
+    python $code -i $TILED -o SpaceTx
+    """
+}
+
+process starfish_decode {
+
+    // Use the container specification from the parameter file
+    // No change to this line is required
+    container "${params.contPfx}${module.container}:${module.version}"
+
+    // Specify the project subdirectory for writing the outputs to
+    // The pattern: specification must match the output: files below
+    // TODO: replace report with the desired output directory
+    // TODO: replace the pattern to match the output: clause below
+    publishDir "${params.in}/iss_processing", mode: 'copy', pattern: "*.csv"
+
+    // Stores .command.sh and .command.log from the work directory
+    //   to the project provenance
+    // No change to this line is required
+    publishDir "${Flow.QC(params.in, 'provenance')}", mode: 'copy',
+      pattern: '.command.{sh,log}',
+      saveAs: {fn -> fn.replace('.command', "${module.name}-${task.index}")}
+
+    // Inputs for the process
+    // mcp - MCMICRO parameters (workflow, options, etc.)
+    // module - module specifications (name, container, options, etc.)
+    // img/sft - pairs of images and their matching spatial feature tables
+  input:
+    val mcp
+    val module
+    path code
+    path SpaceTx
+
+    // Process outputs that should be captured and
+    //  a) returned as results
+    //  b) published to the project directory
+    // TODO: replace *.html with the pattern of the tool output files
+  output:
+    path("*.csv"), emit: results
+
+    // Provenance files -- no change is needed here
+    tuple path('.command.sh'), path('.command.log')
+
+    // Specifies whether to run the process
+    // Here, we simply take the flag from the workflow parameters
+    // TODO: change snr to match the true/false workflow parameter in defaults.yml
+  when: mcp.workflow["iss_decoding"]
+
+    // The command to be executed inside the tool container
+    // The command must write all outputs to the current working directory (.)
+    // Opts.moduleOpts() will identify and return the appropriate module options
+    """
+    python $code -i SpaceTx
+    """
+}
+
+workflow starfish {
+
+    // Inputs:
+    // mcp - MCMICRO parameters (workflow, options, etc.)
+    // imgs - images
+    // cbk - Codebook
+  take:
+    mcp
+    img_dir
+
+  main:
+
+    // Apply the process to each (image, sft) pair
+
+    // code_tile = Channel.fromPath("$projectDir/starfish/bin/decoding.py")
+    ///code_convert = Channel.fromPath("$projectDir/starfish/bin/decoding.py")
+
+    code_tile = Channel.fromPath("$projectDir/starfish/bin/tiling.py")
+    TILED = starfish_tile(mcp, mcp.modules['iss_decoding'], code_tile, img_dir)
+
+    println TILED
+    code_convert = Channel.fromPath("$projectDir/starfish/bin/format_to_spacetx.py")
+    SpaceTx = starfish_convert(mcp, mcp.modules['iss_decoding'], code_convert, TILED[0])
+
+    code_decode = Channel.fromPath("$projectDir/starfish/bin/decoding.py")
+    results = starfish_decode(mcp, mcp.modules['iss_decoding'], code_decode, SpaceTx[0])
+    // Return the outputs produced by the tool
+
+  emit:
+    results[0]
+}

starfish/bin/decoding.py

@@ -0,0 +1,93 @@
+import os
+import starfish
+import argparse
+from starfish.image import ApplyTransform, LearnTransform, Filter
+from starfish.types import Axes
+from starfish import data, FieldOfView
+from starfish.spots import FindSpots, DecodeSpots, AssignTargets
+
+codebook = data.ISS(use_test_data=True).codebook
+
+
+def get_args():
+    """Get command-line arguments"""
+
+    parser = argparse.ArgumentParser(
+        description='Input/output directories for data formatting',
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+
+    parser.add_argument('-i',
+                        '--input_dir',
+                        default='Tiled',
+                        type=str,
+                        help='Input root directory')
+
+    return parser.parse_args()
+
+
+def iss_pipeline(fov, codebook):
+    #fov = experiment.fov()
+    primary_image = fov.get_image(FieldOfView.PRIMARY_IMAGES)
+    anchor = primary_image.sel({Axes.ROUND: 0})
+    anchor_dots = anchor.reduce({Axes.CH, Axes.ZPLANE}, func="max")
+
+    learn_translation = LearnTransform.Translation(reference_stack=anchor_dots,
+                                                   axes=Axes.ROUND, upsampling=100)
+
+    transforms_list = learn_translation.run(
+        primary_image.reduce({Axes.CH, Axes.ZPLANE}, func="max"))
+
+    warp = ApplyTransform.Warp()
+    registered = warp.run(primary_image, transforms_list=transforms_list, in_place=False, verbose=True)
+
+    # Filter raw data
+    masking_radius = 15
+    filt = Filter.WhiteTophat(masking_radius, is_volume=False)
+    filtered = filt.run(registered, verbose=True, in_place=False)
+    print(filtered)
+
+    bd = FindSpots.BlobDetector(
+        min_sigma=1,
+        max_sigma=3,
+        num_sigma=30,
+        threshold=0.01,
+        measurement_type='mean'
+    )
+
+    # Check if experiment has anchor or not
+    # Locate spots in a reference image:
+    # Old one: spots = bd.run(reference_image=fov.get_image("anchor_dots"), image_stack=filtered)
+    spots = bd.run(reference_image=anchor_dots, image_stack=filtered)
+
+    # decode the pixel traces using the codebook
+    decoder = DecodeSpots.PerRoundMaxChannel(codebook=codebook)
+    decoded = decoder.run(spots=spots)
+
+    return decoded
+#imshow_plane(dots)
+
+# process all the fields of view:
+def process_experiment(experiment: starfish.Experiment, cb: starfish.Codebook):
+    decoded_intensities = {}
+    for i, (name_, fov) in enumerate(experiment.items()):
+        print(name_)
+        decoded = iss_pipeline(fov, codebook=cb)
+        decoded_intensities[name_] = decoded
+
+    return decoded_intensities
+
+def save_fov(fov, dataframe):
+    fov_df = dataframe.to_features_dataframe()
+    fov_df.to_csv(f'fov_{fov}.csv', index=False)
+
+
+def main():
+    args = get_args()
+    exp = starfish.Experiment.from_json(os.path.join(args.input_dir, 'primary', 'experiment.json'))
+    results = process_experiment(exp, codebook)
+    for fov, data in results.items():
+        save_fov(fov, data)
+
+
+if __name__ == '__main__':
+    main()