Wiki
Clone wikianima / Applications HCS Segmentation
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) }
Updated