diff --git a/CHANGES.md b/CHANGES.md index 9745086e..29fbd1d3 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,35 @@ + +### 2024-08-26 + +* Add staging step before illumination - `start-at: staging` +* `phenoimager2mc` added as a staging module - `staging-method: phenoimager2mc` +* illumination is run by default if `start-at: staging` is given and can be toggled off with `illumination: false` +* Add max projection option to recyze on multiple nuclear and/or membrane channels provided - output will have nuclear channel 0 and membrane channel 1 +* `segmentation-nuclear-channel` and `segmentation-membrane-channel` should only be used if segmentation-max-projection is toggled, otherwise the output of recyze would only contain the channels provided in `segmentation-channel` + * `segmentation-max-projection: true` + * `segmentation-channel: 1 3 8 9` + * `segmentation-nuclear-channel: 1 8` + * `segmentation-membrane-channel: 3 9` + +Example parameter file for running staging, registration, segmentation with cellpose using the above options and quantification: +``` +workflow: + start-at: staging + stop-at: quantification + illumination: false + staging-method: phenoimager2mc + segmentation-recyze: true + segmentation-max-projection: true + segmentation-channel: 1 5 7 8 11 + segmentation-nuclear-channel: 5 11 + segmentation-membrane-channel: 1 7 8 + segmentation: cellpose +options: + phenoimager2mc: -m 6 --normalization max + ashlar: --align-channel 4 --flip-y + cellpose: --pretrained_model cyto --chan 1 --chan2 0 --no_npy +``` + ### 2024-03-10 * Allow for dynamic sample name specification in the fileseries/filepattern expressions. diff --git a/config/defaults.yml b/config/defaults.yml index 188adf7a..70fe72d2 100644 --- a/config/defaults.yml +++ b/config/defaults.yml @@ -4,12 +4,15 @@ workflow: qc-files: inherit tma: false viz: false + illumination: true background: false background-method: backsub + staging-method: phenoimager2mc multi-formats: '{.xdce,.nd,.scan,.htd}' single-formats: '{.ome.tiff,.ome.tif,.rcpnl,.btf,.nd2,.tif,.czi}' segmentation: unmicst segmentation-recyze: false + segmentation-max-projection: false downstream: scimap options: ashlar: -m 30 @@ -33,6 +36,15 @@ modules: version: 2.2.9 cmd: python /app/UNetCoreograph.py --outputPath . input: --imagePath + staging: + - + name: phenoimager2mc + container: ghcr.io/schapirolabor/phenoimager2mc + version: v0.1.1 + cmd: |- + python /phenoimager2mc/scripts/phenoimager2mc.py \ + --indir ${indir} \ + -o ${cycle}.ome.tif background: - name: backsub diff --git a/config/schema.yml b/config/schema.yml index c0723b53..7d905d2d 100644 --- a/config/schema.yml +++ b/config/schema.yml @@ -7,13 +7,18 @@ workflow: - viz - multi-formats - single-formats + - illumination - segmentation - segmentation-channel - segmentation-recyze + - segmentation-nuclear-channel + - segmentation-membrane-channel + - segmentation-max-projection - downstream - ilastik-model - mesmer-model - background + - staging-method - cellpose-model - background-method - flowsom-model diff --git a/docs/parameters/core.md b/docs/parameters/core.md index 03a0eea9..b85db96b 100755 --- a/docs/parameters/core.md +++ b/docs/parameters/core.md @@ -111,6 +111,7 @@ By default, MCMICRO skips this step as it requires manual inspection of the outp ``` yaml workflow: start-at: illumination + illumination: true ``` * Running outside of MCMICRO: [Instructions](https://github.com/labsyspharm/basic-illumination#running-as-a-docker-container){:target="_blank"}. diff --git a/docs/parameters/other.md b/docs/parameters/other.md index 9cce2bb5..9b95122c 100755 --- a/docs/parameters/other.md +++ b/docs/parameters/other.md @@ -7,6 +7,8 @@ nav_exclude: true --- # Other Modules +Staging +1. [phenoimager2mc](./other.html#phenoimager2mc) Segmentation 1. [Ilastik](./other.html#ilastik) @@ -27,6 +29,39 @@ Background subtraction --- +## Staging - Phenoimager2mc +{: .fw-500} + +### Description +Introducing an additional `staging` step in the pipeline, the [phenoimager2mc](https://github.com/schapirolabor/phenoimager2mc){:target="_blank"} module takes in individual unmixed component data tiles produced by the [InForm software by Akoya](https://www.akoyabio.com/phenoimager/inform-tissue-finder/) and produces an `ome-tif` file per cycle that is compatible with ASHLAR. + +### Usage +By default, `staging` is not performed and the parameter has to be provided as shown below. In addition, the `staging-method: phenoimager2mc` parameter can specify which staging option should be used. It should be noted, `illumination` is run by default after `staging` and should actively be turned off if not needed as presented. It is highly recommended that the input tiles already have overlaps between them, if not, gaps will be introduced. + +* Example `params.yml`: + +``` yaml +workflow: + start-at: staging + staging: true + illumination: false + staging-method: phenoimager2mc +options: + phenoimager2mc: -m 6 --normalization max +``` +* Specify number of channels per cycle: `-m` +* Specify normalization (float32 -> uint16) method: `--normalization`, options `max`, `99th` for maximum value normalization or 99th percentile, respectively. +* Running outside of MCMICRO: [Instructions](https://github.com/labsyspharm/mcmicro-ilastik){:target="_blank"}. + + +### Input +Inputs should be saved in the `staging` directory as one subdirectory per cycle containing the `tif` or `tiff` component data files. The `marker` file should represent the expected registered image, file paths should not contain whitespaces, and it should be taken into account that cycle order is the alphabetical order of subdirectory names. + +### Output +An normalized `.ome.tif` file compatible with ASHLAR containing information from all tiles per cycle. + +[Back to top](./other.html#other-modules){: .btn .btn-purple} + ## Ilastik {: .fw-500} @@ -158,16 +193,31 @@ Cellpose is a DL segmentation algorithm able to segment the nuclear or cytoplasm To use this segmentation method add the line `segmentation: cellpose` in the workflow section of the `params.yml` file. Under the options section of `params.yml` specify the input arguments of the cellpose script, such as segmentation model and channel(s) on which the model will be applied. Notice that the channel(s) argument(s), i.e. --chan and --chan2, expect a zero-based index. -For large data sets it is recommended to use the parameters `segmentation-recyze: true` along with `segmentation-channel:`. In the example below we consider an image stack of 10 channels with the nuclear marker in channel 2 and membrane marker in channel 7. The use of `segmentation-recyze: true` will reduce the image stack to these two channels prior to segmentation, hence reindexing the stack channels such that 2-->0 and 7-->1. +For large data sets it is recommended to use the parameters `segmentation-recyze: true` along with `segmentation-channel:`. In the example below we consider an image stack of 12 channels with the nuclear marker in channel 2 and membrane marker in channel 7. The use of `segmentation-recyze: true` will reduce the image stack to these two channels prior to segmentation, hence reindexing the stack channels such that 2-->0 and 7-->1. +If `segmentation-max-projection: true` and multiple channels are provided for nuclear and membrane stains with `segmentation-nuclear-channel:` and `segmentation-membrane-channel:`, the returned image stack will contain the maximum projection of the respective nuclear (recyze output channel 0) and membrane channels (recyze output channel 1). * Example `params.yml`: ``` yaml workflow: - segmentation-channel: 2 7 + segmentation-channel: 2 7 + segmentation-recyze: true + segmentation: cellpose +options: + cellpose: --pretrained_model cyto --chan 1 --chan2 0 --no_npy +``` + +* Example `params.yml` using max-projection - the input image is large, so `segmentation-recyze` is toggled. The segmentation in this case should be run on the maximum projection of two nuclear markers (channels 5 and 11), and 3 membrane markers (channels 2 3 and 7). The output from `Recyze` will be a two-channel `tif`, the nuclear channel max projection - channel 0, is provided to the `cyto` model with `--chan2 0`, and the membrane channel max projection - channel 1, is provided to the model with `--chan 1`. + +``` yaml +workflow: + segmentation-channel: 2 3 5 7 11 segmentation-recyze: true segmentation: cellpose + segmentation-max-projection: true + segmentation-nuclear-channel: 5 11 + segmentation-membrane-channel: 2 3 7 options: cellpose: --pretrained_model cyto --chan 1 --chan2 0 --no_npy ``` diff --git a/lib/mcmicro/Flow.groovy b/lib/mcmicro/Flow.groovy index 776eb040..2c009d77 100644 --- a/lib/mcmicro/Flow.groovy +++ b/lib/mcmicro/Flow.groovy @@ -26,15 +26,16 @@ static def flowSegment(wfp) { // Valid start/stop steps in the mcmicro pipeline List mcsteps = [ - "raw", // Step 0 - "illumination", // Step 1 - "registration", // Step 2 - "background", // Step 3 - "dearray", // Step 4 - "segmentation", // Step 5 - "watershed", // Step 6 - "quantification", // Step 7 - "downstream" // Step 8 + "staging", // Step 0 + "raw", // Step 1 + "illumination", // Step 2 + "registration", // Step 3 + "background", // Step 4 + "dearray", // Step 5 + "segmentation", // Step 6 + "watershed", // Step 7 + "quantification", // Step 8 + "downstream" // Step 9 ] // Identify starting and stopping indices @@ -46,7 +47,7 @@ static def flowSegment(wfp) { throw new Exception("Unknown stopping step ${wfp['stop-at']}") // Advance segmentation -> watershed to ensure no dangling probability maps - if( idxStop == 5 ) idxStop = 6 + if( idxStop == 6 ) idxStop = 7 return [idxStart, idxStop] } @@ -62,14 +63,15 @@ static def precomputed(wfp) { // Define whether a precomputed intermediate is relevant [ - raw: idxStart <= 2, - illumination: idxStart == 2, - registration: idxStart == 3 || (idxStart == 4 && !wfp.background) || (idxStart > 4 && !wfp.tma && !wfp.background), // needed for background (3), tma if no background (4), everything else if both tma and background aren't specified + staging: idxStart == 0, + raw: idxStart <= 3 && idxStart > 0, + illumination: idxStart == 3, + registration: idxStart == 4 || (idxStart == 5 && !wfp.background) || (idxStart > 5 && !wfp.tma && !wfp.background), // needed for background (3), tma if no background (4), everything else if both tma and background aren't specified background: idxStart > 3 && wfp.background, // if background specified, required - dearray: idxStart > 4 && wfp.tma, // if tma specified, required - 'probability-maps': idxStart == 6, - segmentation: idxStart == 7, - quantification: idxStart == 8 + dearray: idxStart > 5 && wfp.tma, // if tma specified, required + 'probability-maps': idxStart == 7, + segmentation: idxStart == 8, + quantification: idxStart == 9 ] } @@ -84,22 +86,24 @@ static def doirun(step, wfp) { def (idxStart, idxStop) = flowSegment(wfp) switch(step) { + case 'staging': + return(idxStart == 0 && idxStop >= 0) case 'illumination': - return(idxStart <= 1 && idxStop >= 1) + return(idxStart <= 2 && idxStop >= 2 && wfp.illumination) case 'registration': - return(idxStart <= 2 && idxStop >= 2) + return(idxStart <= 3 && idxStop >= 3) case 'background': - return(idxStart <= 3 && idxStop >= 3 && wfp.background) + return(idxStart <= 4 && idxStop >= 4 && wfp.background) case 'dearray': - return(idxStart <= 4 && idxStop >= 4 && wfp.tma) + return(idxStart <= 5 && idxStop >= 5 && wfp.tma) case 'segmentation': - return(idxStart <= 5 && idxStop >= 5) - case 'watershed': return(idxStart <= 6 && idxStop >= 6) - case 'quantification': + case 'watershed': return(idxStart <= 7 && idxStop >= 7) - case 'downstream': + case 'quantification': return(idxStart <= 8 && idxStop >= 8) + case 'downstream': + return(idxStart <= 9 && idxStop >= 9) case 'viz': return(wfp.viz) default: diff --git a/lib/mcmicro/Opts.groovy b/lib/mcmicro/Opts.groovy index 7d67414f..d1299e9f 100644 --- a/lib/mcmicro/Opts.groovy +++ b/lib/mcmicro/Opts.groovy @@ -138,6 +138,24 @@ static def validateWFParams(wfp, fns) { "segmentation-channel provided" throw new Exception(msg) } + if(!wfp['segmentation-max-projection'] && + (wfp.containsKey('segmentation-nuclear-channel') || wfp.containsKey('segmentation-membrane-channel'))) { + String msg = "Multiple nuclear or membrane channels were requested " + + "but no maximum projection specification is provided. " + + "Either add the segmentation-max-projection parameter " + + "or only use segmentation-channel for channel selection." + throw new Exception(msg) + } + if(wfp['segmentation-max-projection'] && + !(wfp.containsKey('segmentation-nuclear-channel') || wfp.containsKey('segmentation-membrane-channel'))) { + String msg = "Maximum projection specification provided but no " + + "nuclear or membrane channels defined. " + + "Either specify multiple nuclear (and membrane channels) with " + + "segmentation-nuclear-channel (and segmentation-membrane-channel) " + + "or exclude segmentation-max-projection and only use segmentation-channel " + + "for channel specification." + throw new Exception(msg) + } } /** @@ -188,6 +206,18 @@ static def parseParams(gp, fns, fnw) { else mcp.modules['background'] = mcp.modules['background'][0] + // Select the staging module based on --staging-method + mcp.modules['staging'] = mcp.modules['staging'].findAll{ + it.name == mcp.workflow['staging-method'] + } + if(mcp.modules['staging'].size() < 1) { + String msg = "Unknown staging method " + + mcp.workflow['staging-method'] + throw new Exception(msg) + } + else + mcp.modules['staging'] = mcp.modules['staging'][0] + // Filter segmentation modules based on --segmentation mcp.modules['segmentation'] = mcp.modules['segmentation'].findAll{ mcp.workflow.segmentation.contains(it.name) @@ -243,10 +273,61 @@ static def moduleOpts(module, mcp) { copts = module.channel + ' ' + idx.join(' ') } + String ncopts = '' + if(wfp.containsKey('segmentation-nuclear-channel') && + module.containsKey('nuclear-channel')) { + + // Module spec must specify whether indexing is 0-based or 1-based + if(!module.containsKey('idxbase')) + error module.name + " spec in modules.yml is missing idxbase key" + + // Identify the list of indices + List idx = wfp['segmentation-nuclear-channel'].toString().tokenize() + + // Account for recyze, if appropriate + if(wfp['segmentation-recyze']) + idx = (1..idx.size()).collect{it} + + // Account for 0-based indexing + if(module.idxbase == 0) + idx = idx.collect{"${(it as int)-1}"} + + // S3segmenter will work with the first index only + if(module.name == 's3seg') + idx = idx[0..0] + + ncopts = module.nuclear_channel + ' ' + idx.join(' ') + } + + String mcopts = '' + if(wfp.containsKey('segmentation-membrane-channel') && + module.containsKey('membrane-channel')) { + + // Module spec must specify whether indexing is 0-based or 1-based + if(!module.containsKey('idxbase')) + error module.name + " spec in modules.yml is missing idxbase key" + + // Identify the list of indices + List idx = wfp['segmentation-membrane-channel'].toString().tokenize() + + // Account for recyze, if appropriate + if(wfp['segmentation-recyze']) + idx = (1..idx.size()).collect{it} + + // Account for 0-based indexing + if(module.idxbase == 0) + idx = idx.collect{"${(it as int)-1}"} + + // S3segmenter will work with the first index only + if(module.name == 's3seg') + idx = idx[0..0] + + mcopts = module.membrane_channel + ' ' + idx.join(' ') + } // Identify all remaining module options String mopts = '' if(mcp.options.containsKey(module.name)) mopts = mcp.options[module.name] - copts + ' ' + mopts + copts + ' ' + ncopts + ' ' + mcopts + ' ' + mopts } diff --git a/lib/mcmicro/Util.groovy b/lib/mcmicro/Util.groovy index 66f32ab3..0e1e9d3b 100644 --- a/lib/mcmicro/Util.groovy +++ b/lib/mcmicro/Util.groovy @@ -6,7 +6,11 @@ static def getSampleName(f, rawdir) { rel.contains('/') ? rel.split('/').head() : rawdir.parent.getName() } - +static def getCycleNameFromDir(f, rawdir) { + // Resolve paths relative to the input project directory + String rel = rawdir.relativize(f).toString() + rel.split('/').head() +} /** * Extracts a file ID as the first token before delim in the filename diff --git a/main.nf b/main.nf index 507f68d9..1ce37d18 100644 --- a/main.nf +++ b/main.nf @@ -48,7 +48,9 @@ findFiles0 = { key, pattern -> pre[key] ? findFiles = { key, pattern, ife -> pre[key] ? Channel.fromPath("${params.in}/$key/$pattern").ifEmpty(ife) : Channel.empty() } - +findDirs = { key, ife -> pre[key] ? + Channel.fromPath("${params.in}/$key/*", type: 'dir').ifEmpty(ife) : Channel.empty() +} // 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 @@ -59,9 +61,15 @@ findFiles = { key, pattern, ife -> pre[key] ? (formatType, formatPattern) = file("${params.in}/raw/**${wfp['multi-formats']}") ? ["multi", wfp['multi-formats']] : ["single", wfp['single-formats']] -rawFiles = findFiles('raw', "**${formatPattern}", - {error "No images found in ${params.in}/raw"}) - + +stagingDirs = findDirs('staging', + {error "No subdirectories found in staging directory"}) +staging_in = stagingDirs + .map{ tuple( + Util.getSampleName(it, file("${params.in}/staging")), + Util.getCycleNameFromDir(it, file("${params.in}/staging")), + formatType == "single" ? it : it.parent + )} // Here we assemble tuples of 1) path to stage for each raw image (might be a // directory) and 2) relative path to the main file for each image. Processes // must input the first as a path and the second as a val to avoid incorrect or @@ -70,6 +78,8 @@ rawFiles = findFiles('raw', "**${formatPattern}", // used when interpolating these paths into script strings, as we are bypassing // the normal way that paths are passed to channels which handles this escaping // automatically. +rawFiles = findFiles('raw', "**${formatPattern}", + {error "No images found in ${params.in}/raw"}) raw = rawFiles .map{ tuple( Util.getSampleName(it, file("${params.in}/raw")), @@ -105,6 +115,7 @@ pre_qty = findFiles('quantification', "*.csv", {error "No quantification tables in ${params.in}/quantification"}) // Import individual modules +include {staging} from "$projectDir/modules/staging" include {illumination} from "$projectDir/modules/illumination" include {registration} from "$projectDir/modules/registration" include {dearray} from "$projectDir/modules/dearray" @@ -116,6 +127,19 @@ include {background} from "$projectDir/modules/background" // Define the primary mcmicro workflow workflow { + staging(mcp, staging_in, chMrk) + //staging.out.view() + staging.out.map{ + sample, cycle, path -> + tuple(sample, cycle, path, path.toString().split('/').last()) + }.toSortedList { a, b -> a[1] <=> b[1] } + .flatMap() + .map{ + sample, cycle, path, name -> + tuple(sample, path, name) + }.set{ sorted_staging } + raw = raw.mix(sorted_staging) + illumination(wfp, mcp.modules['illumination'], raw) registration(mcp, raw, illumination.out.ffp.mix( pre_ffp ), diff --git a/modules/segmentation.nf b/modules/segmentation.nf index 92bf6708..7e498755 100644 --- a/modules/segmentation.nf +++ b/modules/segmentation.nf @@ -75,7 +75,22 @@ workflow segmentation { chan = mcp.workflow.containsKey('segmentation-channel') ? mcp.workflow['segmentation-channel'].toString() .tokenize().collect{"${(it as int)-1}"}.join(' ') : '0' - recyzeOut = roadie('recyze', recyzeIn.toCut, "--channels $chan", false, '', '' ) + + nuc_chan = mcp.workflow.containsKey('segmentation-nuclear-channel') ? + mcp.workflow['segmentation-nuclear-channel'].toString() + .tokenize().collect{"${(it as int)-1}"}.join(' ') : '0' + + mem_chan = mcp.workflow.containsKey('segmentation-membrane-channel') ? + mcp.workflow['segmentation-membrane-channel'].toString() + .tokenize().collect{"${(it as int)-1}"}.join(' ') : '' + + recyzeOut = roadie('recyze', + recyzeIn.toCut, + ["--channels $chan", + mcp.workflow.containsKey('segmentation-nuclear-channel') ? "--nuclear_channels $nuc_chan" : "", + mcp.workflow.containsKey('segmentation-membrane-channel') ? "--membrane_channels $mem_chan" : "", + mcp.workflow['segmentation-max-projection'] ? "--max_projection" : ""].join(" ").trim(), + false, '', '' ) // Determine IDs of images id_cut = recyzeOut.map{ f -> tuple(Util.getFileID(f, '_crop.ome'), f) } diff --git a/modules/staging.nf b/modules/staging.nf new file mode 100644 index 00000000..9dc6c3cd --- /dev/null +++ b/modules/staging.nf @@ -0,0 +1,58 @@ +import mcmicro.* + +import groovy.text.GStringTemplateEngine + +process phenoimager2mc { + // Use the container specification from the parameter file + container "${params.contPfx}${module.container}:${module.version}" + publishDir "${params.in}/raw", mode: "${params.publish_dir_mode}", + pattern: '*.tif' + + // Provenance + publishDir "${Flow.QC(params.in, 'provenance')}", mode: 'copy', + pattern: '.command.{sh,log}', + saveAs: {fn -> fn.replace('.command', "${module.name}")} + + input: + val mcp + val module + tuple val(samplename), val(cycle), path(indir), path(marker) + + output: + tuple val(samplename), val(cycle), path("*.tif"), emit: img + tuple path('.command.sh'), path('.command.log') + + when: Flow.doirun('staging', mcp.workflow) + + script: + def formatMap = [ + marker: marker, + cycle: cycle, + indir: indir, + ] + + def command = new GStringTemplateEngine() + .createTemplate(module.cmd) + .make(formatMap) + .toString() + + """ + ${command} ${Opts.moduleOpts(module, mcp)} + """ +} + +workflow staging { + take: + mcp // MCMICRO parameters as read by Opts.parseParams() + indir // input directories + marker // marker file + + main: + inputs = indir + .combine(marker) + + phenoimager2mc(mcp, mcp.modules['staging'], inputs) + + emit: + phenoimager2mc.out.img +} diff --git a/nextflow_schema.json b/nextflow_schema.json index a5798644..b8710c69 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -35,6 +35,7 @@ "type": "string", "description": "Name of the first step to be executed by the pipeline.", "enum": [ + "staging", "illumination", "registration", "background", @@ -51,6 +52,7 @@ "description": "Name of the final step to be executed by the pipeline.", "default": "quantification", "enum": [ + "staging", "illumination", "registration", "background", @@ -69,6 +71,14 @@ "type": "boolean", "description": "If selected, MCMICRO treats input data as if background channels should be subtracted using the `Backsub` module. See mcmicro.org documentation for details on additional `background` and `remove` columns required in the `markers.csv` file." }, + "staging": { + "type": "boolean", + "description": "Is a staging step required before stitching with ASHLAR?" + }, + "illumination": { + "type": "boolean", + "description": "Should illumination correction be performed?" + }, "segmentation": { "type": "string", "description": "A list of segmentation modules to run", diff --git a/roadie/scripts/recyze.py b/roadie/scripts/recyze.py index 25d6208f..275ccc20 100644 --- a/roadie/scripts/recyze.py +++ b/roadie/scripts/recyze.py @@ -12,7 +12,7 @@ class PyramidWriter: def __init__( - self, _in_path, _out_path, _channels, _x, _y, _x2, _y2, _w, _h, scale=2, tile_size=1024, peak_size=1024, + self, _in_path, _out_path, _channels, _nuclear_channels, _membrane_channels, _max_projection, _x, _y, _x2, _y2, _w, _h, scale=2, tile_size=1024, peak_size=1024, verbose=False ): if tile_size % 16 != 0: @@ -22,6 +22,7 @@ def __init__( self.in_data = zarr.open(self.in_tiff.series[0].aszarr()) self.out_path = Path(_out_path) self.metadata = from_tiff(self.in_path) + self.max_projection = _max_projection self.tile_size = tile_size self.peak_size = peak_size @@ -29,18 +30,48 @@ def __init__( if self.in_data[0].ndim == 3: # Multi-channel image self.single_channel = False if _channels: - if max(_channels) > self.in_data[0].shape[0]: + if max(_channels) >= self.in_data[0].shape[0]: print("Channel out of range", file=sys.stderr) sys.exit(1) else: self.channels = _channels else: self.channels = np.arange(self.in_data[0].shape[0], dtype=int).tolist() + if _nuclear_channels: + if max(_nuclear_channels) >= self.in_data[0].shape[0]: + print("Nuclear channel out of range", file=sys.stderr) + sys.exit(1) + else: + self.nuclear_channels = _nuclear_channels + else: + self.nuclear_channels = [] + if _membrane_channels: + if max(_membrane_channels) >= self.in_data[0].shape[0]: + print("Membrane channel out of range", file=sys.stderr) + sys.exit(1) + else: + self.membrane_channels = _membrane_channels + else: + self.membrane_channels = [] else: # Single Channel image self.single_channel = True if _channels and max(_channels) > 0: print("Channel out of range", file=sys.stderr) sys.exit(1) + if _nuclear_channels and max(_nuclear_channels) > 0: + print("Nuclear channel out of range", file=sys.stderr) + sys.exit(1) + if _membrane_channels and max(_membrane_channels) > 0: + print("Membrane channel out of range", file=sys.stderr) + sys.exit(1) + if _max_projection: + print("Max projection not possible on single channel image.", file=sys.stderr) + sys.exit(1) + self.channels = [0] + + if self.max_projection and len(self.membrane_channels) > 0 and len(self.nuclear_channels) > 0: + self.channels = [0, 1] + elif self.max_projection and len(self.nuclear_channels) > 0: self.channels = [0] xy = _x is not None and _y is not None @@ -75,7 +106,9 @@ def __init__( print('Params:', 'x', self.x, 'y', self.y, 'height', self.height, 'width', self.width, 'levels', self.num_levels, - 'channels', self.channels) + 'channels', self.channels, 'nuclear_channels', self.nuclear_channels, 'membrane_channels', self.membrane_channels, + 'max_projection', self.max_projection + ) self.verbose = verbose @@ -114,25 +147,39 @@ def tile_shapes(self): return tile_shapes + def max_projection_channel(self, channels): + "Compute the maximum projection of the specified channels." + if not channels: + channels = self.channels + maxprojection = np.max( + self.in_data[0] + .get_orthogonal_selection(( + channels, + slice(self.y,self.y + self.height), + slice(self.x,self.x + self.width ))),axis=0) + return maxprojection + def base_tiles(self): h, w = self.base_shape th, tw = self.tile_shapes[0] - for ci in self.channels: - if self.verbose: - print(f" Channel {ci}:") - if self.single_channel: - img = self.in_data[0][self.y:self.y + self.height, self.x:self.x + self.width] + if self.max_projection: + # If max projection is enabled, generate the maximum projection image + if self.nuclear_channels: + nuclear_img = self.max_projection_channel(self.nuclear_channels) + imgs = [nuclear_img] + if self.membrane_channels: + membrane_img = self.max_projection_channel(self.membrane_channels) + imgs.append(membrane_img) else: - img = self.in_data[0][ci, self.y:self.y + self.height, self.x:self.x + self.width] - print('Shape', img.shape) + imgs = [self.max_projection_channel(self.channels)] + else: + imgs = [self.in_data[0][ci, self.y:self.y + self.height, self.x:self.x + self.width] for ci in self.channels] + + for img in imgs: for y in range(0, h, th): for x in range(0, w, tw): - # Returning a copy makes the array contiguous, avoiding - # a severely unoptimized code path in ndarray.tofile. yield img[y:y + th, x:x + tw].copy() - # Allow img to be freed immediately to avoid keeping it in - # memory while the next loop iteration calls assemble_channel. img = None def cropped_subres_image(self, base_img, level): @@ -241,17 +288,27 @@ def run(self): '--channels', type=int, nargs="+", required=False, default=None, metavar="C", help="Channels to keep (Default: all)", ) + parser.add_argument( + '--nuclear_channels', type=int, nargs="+", required=False, default=None, metavar="N", + help="Specifying nuclear channels to keep", + ) + parser.add_argument( + '--membrane_channels', type=int, nargs="+", required=False, default=None, metavar="M", + help="Specifying membrane channels to keep", + ) + parser.add_argument( + '--max_projection', action='store_true', help="Use max projection", + ) parser.add_argument( '--num-threads', type=int, required=False, default=0, metavar="N", help="Worker thread count (Default: auto-scale based on number of available CPUs)", ) - argument = parser.parse_args() + args = parser.parse_args() # Automatically infer the output filename, if not specified - in_path = vars(argument)['in'] - out_path = argument.out + in_path = vars(args)['in'] + out_path = args.out if out_path is None: - # Tokenize the input filename and insert "_crop" # at the appropriate location tokens = os.path.basename(in_path).split(os.extsep) @@ -264,7 +321,7 @@ def run(self): stem = os.extsep.join(tokens[0:-1]) + "_crop" out_path = os.extsep.join([stem, tokens[-1]]) - num_threads = argument.num_threads + num_threads = args.num_threads if num_threads == 0: if hasattr(os, "sched_getaffinity"): num_threads = len(os.sched_getaffinity(0)) @@ -273,6 +330,5 @@ def run(self): tifffile.TIFF.MAXWORKERS = num_threads tifffile.TIFF.MAXIOWORKERS = num_threads * 5 - writer = PyramidWriter(in_path, out_path, argument.channels, argument.x, argument.y, - argument.x2, argument.y2, argument.w, argument.h) + writer = PyramidWriter(in_path, out_path, args.channels, args.nuclear_channels, args.membrane_channels, args.max_projection, args.x, args.y, args.x2, args.y2, args.w, args.h) writer.run() diff --git a/roadie/tasks.yml b/roadie/tasks.yml index d32bdb25..fd61268a 100644 --- a/roadie/tasks.yml +++ b/roadie/tasks.yml @@ -1,7 +1,7 @@ recyze: description: channel and pixel cropping output: '*.ome.tif*' - params: ['x', 'x2', 'y', 'y2', 'w', 'h', 'channels'] + params: ['x', 'x2', 'y', 'y2', 'w', 'h', 'channels', 'nuclear_channels', 'membrane_channels', 'max_projection'] pyramidize: description: generate a pyramid for the input image output: '*.ome.tif'