Wiki

This pipeline example segments a large set of generated DAPI staining like images. The set contains images with various out-of-focus levels, and numbers of cells.
The pipeline tries three different segmentation methods and parallelizes the segmentation. Each method is compared to the known number of cells in every image.
The pipeline downloads the data directly from Broad institute. Note, that MATLAB is required by the pipeline.
#!/usr/bin/env anduril
//$OPT --pipe "tee $ANDURIL_EXECUTION_DIR/_log"
//$OPT --pipe "anduril-pager --ls -t --title $ANDURIL_EXECUTION_DIR"
//$OPT --threads 12

import anduril.builtin._
import anduril.tools._
import anduril.anima._
import org.anduril.runtime._

object Analysis {

/* This first part defines variables, functions and scripts to be used
 * later in the script
*/
// Download the data set
var filedownload=URLInput(url="https://www.broadinstitute.org/bbbc/BBBC005/BBBC005_v1_images.zip")
// Custom script to extract the archive
val fileunarchive=QuickBash(filedownload, script="mkdir $out; cd $out; unzip -j $in")
val fileselection=QuickBash(in=fileunarchive, script="mkdir $out; cd $out; mv $in/*w1.TIF .")

// Define image input
val filedir=fileselection.out

// Defines how many parts the data is split in to
var parals=24

// Fiji script to watershed a mask
//   Returns the mask in output port 'out'
//   and the perimeter of the mask in port 'out2'
val watershed_script="""

IJ.run(in1, "8-bit", "");
IJ.setThreshold(in1,127, 255)
IJ.run(in1, "Convert to Mask", "")

IJ.run(in1, "Invert", "")
IJ.run(in1, "Fill Holes", "")
IJ.run(in1, "Watershed", "")

out2=imageCopy(in1)
IJ.run(out2, "Outline", "")
IJ.run(out2, "Invert", "")

out=in1
IJ.run(out, "Invert", "")
"""

val watershed_scriptWL="""

IJ.run(in1, "8-bit", "");
IJ.setThreshold(in1,127, 255)
IJ.run(in1, "Convert to Mask", "")

IJ.run(in1, "Fill Holes", "")
IJ.run(in1, "Watershed", "")

out=in1
out2=imageCopy(in1)

IJ.run(out2, "Outline", "")
IJ.run(out2, "Invert", "")
IJ.run(out, "Invert", "")
"""


// This function calculates a difference between the number of 
// segmented objects and actual number
def validation(meas: CSV,summary: Boolean=false): CSV = {
    val measurements_summary=if (!summary) {
        CSVSummary(meas,clusterCol="File").out
    } else {
        meas
    }
    val measurement_validation=CSVTransformer(csv1=measurements_summary,
        transform1="as.integer(gsub('SIMCEPImages_.*_C([0-9]+)_F.*','\\\\1',csv1[,'File']))",
        transform2="csv1[,c('Count')]",
        transform3="as.integer(abs(transformed[,1] - csv1[,'Count']))",
        columnNames="c('RealCells','Count','Error')")
    val validation_summary=CSVSummary(measurement_validation, summaryType="sum", counts=false)
    // Here is the calculation of non-incorrect ratio
    val validation_ratio=CSVTransformer(csv1=validation_summary, transform1="csv1",
                     transform2="round(1000*(1-(csv1[,3]/csv1[,1])))/10", 
                     columnNames="c(colnames(csv1),'Correct ratio')")
    validation_ratio
}
/* In this part, we run the analysis, and parallelize with for-loops
*/

val visuals=NamedMap[ImageList]("visuals")
val numericFJ=NamedMap[CSV]("numericFJ")
val numericFS=NamedMap[CSV]("numericFS")
val numericWL=NamedMap[CSV]("numericWL")
val file_split=FolderSplit(filedir,N=parals,
                              link=true,order="random")

for ( (x,f) <- iterArray(file_split.out) ) {
    // fetch [parals]:th part of the data set (gray scale images)
    val file_set=file_split.out(x).force()

    // Segment the images with Otsu algorithm, and correct the level by a constant.
    // The correction is set by visual evaluation by running only one 
    // parallel thread first.
    val FJnuclmask=SegmentFiji(file_set,
                        method="Otsu",
                        corr=1.3)
    // Watershed the mask
    val FJwaterprocess=Fiji(
                        in1=FJnuclmask.mask,
                        script=watershed_script)
    // Measure the features for each object
    val FJfeatures=FijiFeatures(mask=FJwaterprocess.out,
                                images=file_set, invert=true)

    // Get masks with FarSight library 
    val FSnuclmask=SegmentGraphCut(file_set,
                     maxScale=12,
                     minScale=10)
    // Make sure mask is valid 8-connected for Fiji
    val FSnuclsep=ImageMagick(in1=FSnuclmask.mask,
                     command="@in1@ -morphology Erode Disk:1 @out@")

    // Measure the features for each object
    val FSfeatures=FijiFeatures(mask=FSnuclsep,
                                     images=file_set)
    // Get masks with Wavelet segmentation library (requires MATLAB)
    val WLnuclmask=SegmentImage(file_set,
                                method="Wavelet", clearBorders=false, fill=false,
                                waveletParam="4,4")
    // Watershed the mask
    val WLwaterprocess=Fiji(
                        in1=WLnuclmask.mask,
                        script=watershed_scriptWL)

    val WLfeatures=FijiFeatures(mask=WLwaterprocess.out,
                                     images=file_set,invert=true)
    // Create an RGB visualization of the grayscale image, and the 
    // segmentation perimeter. Change to PNG for smaller file size.
    // For brevity, these are done with ImageMagick, instead of Fiji
    var compose_command="-compose Plus @in1@ \\( @in2@ +level-colors ,Green \\) -composite @out@"
    val FJvisu=ImageMagick(in1=file_set, in2=FJwaterprocess.out2,
                    oldExtension=".TIF", extension=".png",
                    command=compose_command)
    val FSvisu=ImageMagick(in1=file_set, in2=FSnuclmask.perimeter,
                    oldExtension=".TIF", extension=".png",
                    command=compose_command)
    val WLvisu=ImageMagick(in1=file_set, in2=WLwaterprocess.out2,
                    oldExtension=".TIF", extension=".png",
                    command=compose_command)
    val visu=ImageMagick(in=makeArray(FJvisu,FSvisu,WLvisu), binary="montage",
                    command="@inall@ -tile x1 -geometry +3+0 @out@")

    // Remove objects where pixel area is less than 200
    val FJfeatures_filt=CSVFilter(FJfeatures.out, lowBound="Area=200")
    val FSfeatures_filt=CSVFilter(FSfeatures.out, lowBound="Area=200")
    val WLfeatures_filt=CSVFilter(WLfeatures.out, lowBound="Area=200")
    // Store results for later joining.
    visuals(x)=visu.out
    numericFJ(x)=FJfeatures_filt.out
    numericFS(x)=FSfeatures_filt.out
    numericWL(x)=WLfeatures_filt.out

}

val measurements=NamedMap[CSV]("measurements")
measurements("FJ")=CSVListJoin(in=numericFJ, fileCol="")
measurements("FS")=CSVListJoin(in=numericFS, fileCol="")
measurements("WL")=CSVListJoin(in=numericWL, fileCol="")

val method_rename=NamedMap[CSV]("method_rename")
for (s<-measurements) {
    // Parse the filename to find what is the simulated off-focus level
    // This is the final result from the analysis    
    val measurements_grouped=CSVTransformer(csv1=s._2, 
         transform1="csv1",
         transform2="as.integer(gsub('SIMCEPImages_.*_F([0-9]+)_s.*','\\\\1',csv1[,'File']))",
         columnNames="c(colnames(csv1),'Focus')")
    // Split the results by off-focus level
    val measurements_split=CSVSplit(measurements_grouped, labelCol="Focus")
    // Validate each of the focus levels separately
    val validations=NamedMap[CSV]("validations"+s._1)
    for (el<-iterArray(measurements_split.out)) {
        validations(el._1)=validation(meas=measurements_split.out(el._1))
    }    
    // Join validations in one table
    val method=CSVListJoin(in=validations, fileCol="Focus")
    // Sort CSV by Focus level
    val method_sorted=CSVSort(method)
    method_rename(s._1)=CSVFilter(method_sorted.out, 
                         rename="Correct ratio="+s._1)
}

// Join results
val methods_compare=CSVJoin(in=method_rename,
                         intersection=false,
                         keyColumnNames="Focus")

// Plot a graph with non-incorrect ratios.
val  focus_plot=GNUPlot(in=makeArray(methods_compare),
             script="""
    set terminal png font 'Helvetica,11' size 1024,768
    set output @output_segmentation_compare.png@
    set title "Cell segmentation result comparison"
    set ylabel "Ratio of non-incorrect cells [%]"
    set xlabel "Simulated out-of-focus level"
    set key right bottom
    set yrange [60:105]
    set style line 1 lc rgb '#000000' lt 1 lw 2 pt 7 ps 0.7
    set style line 2 lc rgb '#ff0000' lt 1 lw 2 pt 4 ps 0.7
    set style line 3 lc rgb '#0000ff' lt 1 lw 2 pt 6 ps 0.7
    plot @input_1@ using 1:5 title "Fiji, Otsu + Watershed" with lp ls 1,\
          @input_1@ using 1:6 title "FARSIGHT, Graph cut" with lp ls 2,\
          @input_1@ using 1:7 title "MATLAB, Wavelet segmentation + Watershed" with lp ls 3""")

// Join the visualization images from the parallel sets
val visualization=Array2Folder(visuals, fileMode=".", link=false)

val highlight_self=SyntaxHighlight(in = INPUT("hcs2.scala"),
                    language        = "Anduril",
                    rowNumbers      = true,
                    title           = "Code listing",
                    PNGFont         = "Andale Mono",
                    PNGSize         = 14,
                    PNGPad          = 5)

}