# Bioimage Analysis Wiki

### Sidebar

EMBL BioImage Data Analysis

EuBIAS

NEUBIAS

Popularity Ranking

 RT @KuglerElisabeth: One question I often get is "How do I assess my segmentation?" Together with Andrik Bin Rampun, .@Timchico, and Paul… About 3 hours, 9 mins ago by: Kota Miura (@cmci_) RT @clairebrown514: @cmci_ We use https://t.co/J9y0625zRl . People still need to be convinced and reminded to put the #RRID in the acknowled… About 4 hours, 9 mins ago by: Kota Miura (@cmci_) RT @Coronin: Wonderful historical perspective from Tom Pollard: Landmarks in the discovery of a role for actin in cell locomotion https://t… About 12 hours, 45 mins ago by: Kota Miura (@cmci_) RT @IAMichaelNelson: @cmci_ @CCI_liv I recommended to my old core https://t.co/mOai5Zez5b . A "how to use our core facility" protocol and ma… About 18 hours, 8 mins ago by: Kota Miura (@cmci_) RT @cmci_: @CCI_liv It’s an idea. What would you say? There are journals that allows versioning, so it can be updated every time when there… About 18 hours, 26 mins ago by: Kota Miura (@cmci_) @CCI_liv It’s an idea. What would you say? There are journals that allows versioning, so it can be updated every ti… https://t.co/xz3xKR243j About 18 hours, 36 mins ago by: Kota Miura (@cmci_)
documents:140326jointcourse

# BioIT-EMBL Centres Joint Course 2014

See course website (internal) for its aim and outline of the contents:

# Analysis of High-Throughput Microscopy Image Data

BioIT-EMBL Centres Joint Training Week, Mar. 26, 2014

Place: ATC Computer Teaching Lab

Kota Miura (CMCI, EMBL)

## Instructors

• Kota Miura

• Centre for Molecular and Cellular Imaging (CMCI), EMBL
• miura@embl.de
• Christoph Schiklenk

• Haering Lab, Cell Biology & Biophysics, EMBL
• Christoph has been using Fiji/Jython for a year and he now is a quite literate Fiji scripting guy. He is volunteering to assist your learning.
• Clemens Lakner

• Bio-IT, EMBL
• Clemens wrote the job-array generator scripts for image analysis. He will explain to you how his shell script is working.
• Christian Tischer

• ALMF, EMBL
• Tischi will be partially appearing to assist you. Imaging expert and not only an Jython-Fiji expert, but also a Cell Profiler expert.
• Aliaksandr Halavatyi

• Pepperkok Lab, Cell Biology & Biophysics, EMBL
• Aliaksandr will be partially appearing to assist you. He has been coding image analysis tools in Matlab, but knows ImageJ/Fiji pretty well and coding in Jython.
• Acknowledgements

• Thanks to Fatima Verissimo, Faba Neumann, Rainer Pepperkok and Jean-Karim HÃ©richÃ© for providing data, sharing their experiences and giving advices on these data.
• Thanks to Bernd Klaus, Matt Rogon and Aidan Budd (the course team) for their effort in starting up this course and working together.

## Course Outline

1. Aim: Learn how to process and analyze high content image data by Fiji scripting
2. Introduction: the background of data, setting up.
3. Processing and analyzing using Fiji GUI
• We first try to segment nucleus image and do measurement using GUI, to know how the algorithm works.
4. Python Crash Course
• We want to write a script that does processing and measurements: for this we first need to learn basics of python syntax.
5. Fiji scripting in Jython
• Jython is a Java-implemented Python, allows to use full functionality of Fiji by scripting. We learn how to write Jython script to do our task automatically.
6. Submitting Jobs (executing Jython script) to the cluster.
• As there are many images, we use cluster to parallelize the computation, to get results in a short time.

## Introduction

Image based screening (high content screening) is one of the modern methods for screening genes involved in a biological process. The main idea is to systematically acquire a huge amount of data with various treatments, and try to isolate treatments that affects the target biological process.

Dataset we use in this course is from a secretory pathway screening. Proteins synthesised in endoplasmic reticulum are transported first to the Golgi, where proteins are added with modifications, and then to the plasma membrane or other destined intracellular structures.

VSVG is a model protein used for studying this transport system and is a heat-sensitive protein. This protein could be accumulated in / released from ER by temperature shifting, to experimentally analyze the process of intracellular transport from ER to the plasma membrane.

To study the transport process, it is desirable to take time lapse sequence to follow the changes in VSVG localization within cell. However, it is sufficient to assess the transport process by taking a snapshot at a defined time point after releasing the protein from ER and check where the protein is located. In control condition, proteins should be localized in the plasma membrane after successful transports. If the transport out from the ER is interfered due to suppression of the component of transport, proteins should be stably located in the ER. Likewise, finding proteins in the Golgi suggests that transport machinery from Golgi to the plasma membrane is blocked.

With this snapshot strategy, it is possible to systematically prepare a set of siRNAs and drug treatments to test their effect on vesicle transport system: the high content image screening.

In this course, we will try to analyze all images from different types of treatments automatically, evaluate and compare the difference in the efficiency of transport. This COULD be done manually, but we have huge number of various treatments: approximately 700,000 of images to analyze. To not to end the rest of our life clicking images, we will completely automate the image processing and analysis to extract transport parameters and output results as a huge list of parameter table.

This table then is studied using statistical techniques to distinguish treatments that are largely affecting the transport system.

Followings are the literatures published with this dataset.

• Liebel, U., Starkuviene, V., Erfle, H., Simpson, J.C., Poustka, A., Wiemann, S., Pepperkok, R., 2003. A microscope-based screening platform for large-scale functional protein analysis in intact cells. FEBS Lett. 554, 394â€“8. http://www.ncbi.nlm.nih.gov/pubmed/17275941

• Simpson, J.C., Cetin, C., Erfle, H., Joggerst, B., Liebel, U., Ellenberg, J., Pepperkok, R., 2007. An RNAi screening platform to identify secretion machinery in mammalian cells. J. Biotechnol. 129, 352â€“65. http://www.ncbi.nlm.nih.gov/pubmed/17275941

• Simpson, J.C., Joggerst, B., Laketa, V., Verissimo, F., Cetin, C., Erfle, H., Bexiga, M.G., Singan, V.R., HÃ©richÃ©, J.-K., Neumann, B., Mateos, A., Blake, J., Bechtel, S., Benes, V., Wiemann, S., Ellenberg, J., Pepperkok, R., 2012. Genome-wide RNAi screening identifies human proteins with a regulatory function in the early secretory pathway. Nat. Cell Biol. 14, 764â€“74. http://www.ncbi.nlm.nih.gov/pubmed/22660414

## Configuration for March 26, 2014

Here is what you should do to use Fiji locally and via network. We learn GUI and scripting using local Fiji (one that you will install in your machine) and the network Fiji (resides in ALMF server and run it headlessly using command line).

In brief,

• For desktop Fiji, you just need to download the installer.
• For the network Fiji, you need to set path both for the network Fiji and Java virtual machine (Java should be 1.6).

Just a note: if you try to run Fiji from sshint.embl.de, it fails because the Java version there is 1.4.

### Install Fiji on your desktop

Access the following page

http://fiji.sc/Downloads


then install "Life-Line Version", Linux 64-bit. Place the Fiji folder on your desktop.

### Configure access to the Fiji in the network

For using Fiji in the network from command line, do the following.

#### Add an alias in your local bash profile.

Open ~/.bash_profile

open ~/.bash_profile


Add the following line in your bash profile (.bash_profile).

alias cluster='ssh -t <YOUR UNIX USERNAME>@submaster '\''bsub -M 5000 -Is bash'\'''


By the way, this will allow you to execute the command below as an alias 'cluster'

ssh -t <YOUR UNIX USERNAME>@submaster \bsub -M 5000 -Is bash


If you have time, try the above command directly and see that it works.

Reload the .bash_profile by

source .bash_profile


to validate your changes.

#### Add several paths to your .bash_profile in the network.

If your are ready with the above new alias, try

cluster


You will then be asked for your password, and then will be connected to one of the cluster nodes in an interactive mode.

vim ~/.bash_profile


If you are not use to vi editor, you could use others such as

nano ~/.bash_profile



Otherwise, it might be able to edit the .bash_profile directly from your course Linux machine.

Add the following lines in your .bash_profile

export PATH=/g/software/bin:$PATH export PATH=/g/almf/software/bin2:$PATH
export PATH=/g/almf/software/bin/java/jdk1.6.0_45/bin/:$PATH export JAVA_HOME=/g/almf/software/bin/java/jdk1.6.0_45  This is just to run Fiji from your command line (the one that I am maintaining) and to use the right Java version. You could delete these lines after the course. Reload the profile source .bash_profile  ### Check if Fiji is working If you are still connected to one of the nodes, do java -version  The output should look like [miura@compute083]~$
java -version
java version "1.6.0_45"
Java(TM) SE Runtime Environment (build 1.6.0_45-b06)
Java HotSpot(TM) 64-Bit Server VM (build 20.45-b01, mixed mode)


If your Java version is OK (same as above), then

fiji --help


The output should be

export PATH=/g/almf/software/bin/java/jdk1.6.0_45/bin:$PATH  Then source .bashrc  or source .bash_profile  ### Data All data are located under /g/data/bio-it_centres_course/data/VSVG/  Plates are numbered such that plateID-Number--year-month-date  For example: 0001-03--2005-08-01  This means that this is the plate ID 1, third experiment (same plate configuration is repeated at least three times each). We analyze plateID from 1 to 159. There are 703 folders so the cluster job array will be composed of 703 jobs. Within each of these folders, images are stored under data/  ### Example Image Example image set for learning image processing and analysis is in /g/data/bio-it_centres_course/data/course  There are three files from single spot in a plate. ### Scripts The script for processing images for a single plate is /g/almf/software/scripts2/measTransportBatch3.py  To process a plate, for example the first plate '0001-03--2005-08-01', the plate folder name is given as an argument to a command like below: fiji --headless --mem=1000m /g/almf/software/scripts2/measTransportBatch3.py '0001-03--2005-08-01'  Pre-screening script is /g/almf/software/scripts2/listFocusedImagesV2.py  A script for generating job array written by Clemens Lakner is /g/data/bio-it_centres_course/data/VSVG/fiji-sub/fiji_script.sh  The code: <> ## Workflow: GUI The protocol below is described in Simpson (2007). ### Background subtraction Background subtraction will only be done with PM and VSVG images. In the final script, the function name is backgroundSubtraction Get minimum and mean intensity of the image. [Analyze > Set measurementsâ€¦] : min and max, mean [Analyze> Measure]  Results will be shown in the Results Table. We use min and mean values. We call them min1, mean1. Then calculate mean of the pixels those with values greater than min1 and less than the mean1. [Analyze > Set measurementsâ€¦] :select min and max, mean, limit to threshold [Image > Adjust > thresholdâ€¦] click "Set" and put Lower Threshold Level: min1 Upper Threshold level: mean1 [Analyze> Measure] results: mean2 Click "Reset" in Threshold GUI.  Subtract mean2 value from the image. [Process > Math > Subtract]: use "mean2".  ### nucleus segmentation Nucleus segmentation. We only do this with Dapi image. Corresponding function in the script is nucleusSegmentation(imp2). [Image > Type > 8bit] [Process > Filters > Gaussian Blur...] : radius=2.0 [Image . Adjust > Auto Local Threshold] : algorithm, Bernsen, other parameters stay same  Pre-filtering by size and circularity [Analyze > Set Measurements] :select area, standard deviation, shape descriptors, integrated density [Analyze > Analyze Particlesâ€¦] :set area 20 - 1000 :set circularity 0.8 - 1.0 :show: Masks :check Display Results, Clear results, Exclude on Edges, Include Holes  Rename the resulted image as impNuc [Image > Rename]  Duplicate the resulted image and then invert LUT [Image > Duplicateâ€¦] [Image > Lookup Tables > Invert LUT]  Dilate segmented nucleus by 15 pixels (to check if nucleus will merge) [Process > Binary > Options] set iterations to 15, Check "Black Background" [Process > Binary > Dilate]  (called imp2) Particle analysis 1 (with imp2) [Analyze > Analyze Particlesâ€¦] :set area 20 - 10000 :set circularity 0.8 - 1.0 :show: None :check Display Results, Clear results, Exclude on Edges, Include Holes  q1, q3, offset calculation: Copy and paste https://gist.github.com/miura/9641262  Example results: 2926.0 3854.0 1392.0 min2 = q1 - offset (= 2462) max2 = q3 + offset (= 5246)  Particle analysis 2 (remove outlier) with imp2, [Analyze > Analyze Particlesâ€¦] :set area min2 - max2 :set circularity 0.8 - 1.0 :show: Mask :check Display Results, Clear results, Exclude on Edges, Include Holes output: the mask. rename it as impDilatedNuc  Filter original nucleus  [Process > Image Calculator â€¦] AND operation for impNuc and impDilatedNuc  ### Check Segmentation Results.  Visual inspection. We terminate analysis if there is no nucleus segmented.  ### Generate nucleus ROIs from segmented nucleus. Store them as a list.  [Edit > Selection > Create Selection] [Edit > Selection > Make Inverse] [Edit > Selection > To ROImanager]  ### Roi modifications  In ROI manager a. Generate All Cell ROIs (enlarge) Activate the ROI above, then [Image > Selection > Enlargeâ€¦] : 15 pixels click "Add". Second ROI appears in the ROI manager list. this will be the "All cell" b. Generate JuxtaNuclear ROIs (enlarged minus nucleus) select both ROIs, and then [More >> XOR]. click "Add". Third ROI appears in the listing. This ail be the zone surrounding nucleus.  ### Measurements Do measurements using modified ROIs (measureROIs).  select corresponding ROI and [Analyze > Measure] a. enlarged: all cell plasma membrane b. ring rois: juxtanuclear signal for VSVG c. all cell VSVG intensity  ### Transport Efficiency Calculate transport efficiency and add the results to results table. All results from single plate are merged and placed in a file, but we do not do this manually. ## Workflow: Python Primer Jython is Java-implemented Python, allowing us to access Java classes in Python syntax. To write Jython scripts in Fiji, we use Script Editor. [File > New > Script]  This command launches the editor and a new script file. Script Editor has its own menu, and one of them is [Language > Python]  To set the language to Python (Jython). The editor has two panels. Upper one is the input field, and the lower one is the output field. Output field has two types, standard output and error output. You could toggle this by clicking "error" or "output". We learn the basics of Python commands interactively. Below is a list of Python functions we go over. print print with comma separated variables functions return multiple returns list initialization range sorted(list) slicing l.append() len for allcellA = [roiEnlarger(r) for r in nucroiA] enumerate extend() if, else, break String concatenation f.endswith() dictionary initialization d['a'] = 3 Accessing files path / file managements os.listdir os.path.isfile os.path.join Python regex re.compile re.search re.group I/O f = open(path, 'rb') sys.argv csv.reader  Codes used in this interactive session assembled in a file and in https://gist.github.com/miura/9687511 ### File access import os root = '/Volumes/data/bio-it_centres_course/data/course' #filename = 'dapi.tif' files = os.listdir(root) print 'file 1: ', files[0] fullpath = os.path.join(root, files[0]) print fullpath print 'file 2: ', files[1] fullpath = os.path.join(root, files[1]) print os.path.isfile(fullpath)  ### Writing CSV import csv f = open('/Users/miura/tmp/test.csv', 'wb') writer = csv.writer(f) writer.writerow(['id', 'mean', 'sd']) writer.writerow([1, 10.5, 3.3]) f.close()  ### Reading CSV import csv filepath = '/Users/miura/tmp/test.csv' f = open(filepath, 'rb') data = csv.reader(f) for row in data: print ', '.join(row)  ### For more information see http://www.jython.org/docs/library/indexprogress.html http://www.jython.org/docs/library/functions.html ## Workflow: Fiji Jython Scripting ### IJ.log Let's first try accessing Java classes. Try the following. IJ.log("Hello World!")  This will print the argument "Hello World!" in the Log window (in ImageJ macro, similar job can be done by print(string) command). We could add another line IJ.log("Hello World!") IJ.log("\\Clear")  Run this code, and compare with the output of IJ.log("\\Clear") IJ.log("Hello World!")  The double back-slash \\  is called escape sequence, and is a command send to the interpreter. Instead of printing the following string in the Log window, it is interpreted as a command. In case of "\Clear", texts in the Log window will be cleared. ### Javadoc Now, let's see a bit more of details about "IJ.log". Click the following link and see the page Javadoc - IJ This is a documentation called "Javadoc" and is a reference for all the functionality that an Java application has. It's VERY IMPORTANT to understand how to use this documentation, so I wll explain in details. ImageJ is a collection of many classes, and you could find all those classes in this Javadoc. The web page has side bar in the left side, and it is separated to the upper and the lower panels. The upper panel has its title "Packages" and the bottle panel is "Classes". As I have already said, ImageJ is a collection of many classes - and they are all listed in the bottom panel. Since there are many classes, they are organized in hierarchical way by "packages". You could imagine each packages as a folder in file system. For example, you could find a package named "ij". It means something like a folder named "ij", and for a package "ij.plugin.filter", its a folder "filter" under a folder "plguin" within a folder "ij" (in fact, class files are located in file system in this way). In the bottom panel, all classes are listed as long as you have not filtered them by clicking one of the packages in the upper panel. Try to look for the class IJ by scrolling down the lower panel, and click the link "IJ". You will then see a new page in the right side: There is a big text saying "Class IJ". IJ from IJ.log is a class, which contains many methods. To see those methods, scroll down the page and you will find a table titled "Method Summary". This is an alphabetical list of all methods that IJ contains. Among them, you could find log(java.lang.String s)  To use a method from IJ class, you just need a dot after the class name and then add the name of the method. In case of the method "log", you need an argument and it should be a String type. In the example we tried above, it was "Hello World!". How about other methods? Try the method beep(). IJ.beep()  When this code is executed, your machine will make a beep sound. In such a way all the methods are usable. ### Loading and showing Image One of methods available in class IJ is openImage(path). Check it's javadoc information (Javadoc). It says: openImage public static ImagePlus openImage(java.lang.String path) Opens the specified file as a tiff, bmp, dicom, fits, pgm, gif or jpeg image and returns an ImagePlus object if successful. Calls HandleExtraFileTypes plugin if the file type is not recognised. Displays a file open dialog if 'path' is null or an empty string. Note that 'path' can also be a URL. Some reader plugins, including the Bio-Formats plugin, display the image and return null. It loads an image from a file in the path, and returns an ImagePlus object (we will study more about object, or instance in below). imp = IJ.openImage('/g/data/bio-it_centres_course/data/course/--W00005--P00001--Z00000--T00000--dapi.tif') imp.show()  imp now points to the loaded image, and is an instance of class ImagePlus. The ImagePlus holds image itself (pixel data) and metadata (such as its size, image name, units and so on). Check ImagePlus Javadoc ImagePlus It has many methods, and one of them is show(), which creates a windows and displays image on your desktop. ### Class ImagePlus We study what we could get from the instance of ImagePlus object (the imp we got by loading image file). • imp.getTitle() • add a line print imp.getTitle() • imp.duplicate() • add two lines imp2 = imp.duplicate() then imp2.show() • imp2 will be a duplicate of the imp image. • imp.deleteRoi() • imp.setRoi() • imp.getProcessor() ### Class ImageProcessor • impstats = ip.getStatistics() • ip.getMax() • ip.setThreshold(min, max, lut) • ip.resetThreshold() • ip.subtract(val) • ip.invertLut() • r = ip.restRoi() • ip.setRoi() ### Class ImageStatistics • impstats.mean ### Constructing an Object So far, I have not explained an important difference in types of methods: There are static and non-static methods. If you see the Javadoc of class IJ, you will see that all methods of this class has "static" as their property. This means that you could use those methods by appending the name of the method directly after name of the class, such as IJ.log() and IJ.beep(). In contrast, this is not possible with non-static methods. For example, the class ImagePlus has many methods such as "getHeight()", which is to get the height of the image. If you try to do print ImagePlus.getHeight()  It will not work and returns an error. This is simply because we (and mainly the computer) do not know which image that we asking for its height. Instead we must use a specific instance of image, such as the one you retrieve by IJ.openImage(). Remember, imp = IJ.openImage('/Users/miura/image.tif')  This "imp" is a specific instance of image because we know that we loaded the image from a specific file. ImagePlus is a class, but never becomes existing until it is "constructed" by loading an image from a file or creating a image from scratch or grabbing a image that is opened on desktop. Classes that became into existing entity is called an "Instance" or "Object". None static methods work only when an instance is there. Coming back to the static methods of class IJ, their functions are those that does need a presence of specific instance, but only does a general job to any instance or without an instance. In case of IJ.log, we know that there is only one Log window so we do not need to specify an instance. IJ.beep or IJ.openImage as well. In case of images, we know that there could be many different instances of images such as a image from channel 1 and another image from channel2. Though they all belong to a same ImagePlus class, they are different instances of ImagePlus. GaussianBlur is a class that does Gaussian blurring. It's methods are all none static but one, and for these non-static methods to work, an instance of GaussianBlur class should be constructed first. It's also that this class should be imported before it is being in use. Use from [package name] import [class name]  to do so. If this class is not explicitly imported and used, there will be an error message: ImportError: cannot import name GaussianBlur # import class explicitly from ij.plugin.filter import GaussianBlur #Constructs an instance of Gaussian Blur object gb = GaussianBlur()  You could see how to construct an instance in the "constructor summary" in the Javadoc. Now the actual object of GaussianBlur, gb, is available, you could blur an image by using a method blurGaussian(ImageProcessor ip, double sigmaX, double sigmaY, double accuracy)  which would look like gb.blurGaussian(imp.getProcessor(), 5.0, 5.0, 0.02)  The code for loading an image and do GaussianBlur filtering would look like this: #Loading image import os root = '/Volumes/data/bio-it_centres_course/data/course' filename = '--W00005--P00001--Z00000--T00000--dapi.tif' fullpath = os.path.join(root, filename) imp = IJ.openImage(fullpath) # Until here, it's all the same. # import class explicitly from ij.plugin.filter import GaussianBlur # Construct an instance of Gaussian Blur object gb = GaussianBlur() gb.blurGaussian(imp.getProcessor(), 5.0, 5.0, 0.02) imp.show()  Similarly, many classes for processing and analysis works only after constructing its instance. Here is another example, using a class ProfilePlot. Before running this code, please open an image on your desktop, place a line ROI. imp = IJ.getImage() pf = ProfilePlot(imp) profile = pf.getProfile() for val in profile: print val  The second line of this code is the constructor that creates an instance of class ProfilePlot.The constructor in this case uses an argument. See the Javadoc. The class ProfilePlot has three different constructors. • ProfilePlot() • ProfilePlot(ImagePlus imp) • ProfilePlot(ImagePlus imp, boolean averageHorizontally) We could use any of these three constructors. The most basic one would be pf = ProfilePlot()  but then we have not specified an image where we get the intensity profile. This is not usable at all. For this reason, we use pf = ProfilePlot(imp)  to specify the image that we want to work on in the constructor. We now combine the above code with the python package csv to write the output to a comma separated value text file. import csv imp = IJ.getImage() pf = ProfilePlot(imp) profile = pf.getProfile() for val in profile: print val f = open('/Users/miura/Desktop/prof.csv', 'wb') writer = csv.writer(f) for index, val in enumerate(profile): writer.writerow([index, val]) f.close()  ## Workflow: Implementing Image Analysis ### Background Subtraction for VSVG and PM images. Now we start implementing the image processing & analysis of image data. We do this step by step - meaning write several lines, execute to test and then if it's successful, we go on to write more. Let's first set the path to the image file and load it. We've done this already, so you should be able to understand this. import os #path information root = '/g/data/bio-it_centres_course/data/course' filepm = '--W00005--P00001--Z00000--T00000--pm-647.tif' filevsvg = '--W00005--P00001--Z00000--T00000--vsvg-sfp.tif' fullpath = os.path.join(root, filepm) #Load an image imp = IJ.openImage(fullpath) imp.show()  We then get the ImageProcessor object of the current ImagePlus object and do the statistics to get minimum and mean intensity of the current image. #Get the ImageProcessor object of current image ip = imp.getProcessor() #Get ImageStatistics object of current ImageProcessor Object impstats = ip.getStatistics() #Get rtatistics results ipmin = impstats.min ipmean = impstats.mean  getStatistics method returns an ImageStatistics object which does statistical calculations and holds the results. All of these results can be directly accessed as field values (field values - this is new in this tutorial). You could also get other statistics results, which are listed as "Field Summary" of its Javadoc: http://rsbweb.nih.gov/ij/developer/api/ij/process/ImageStatistics.html You could test to print out some of these values. print 'Median:', impstats.median print 'SD:', impstats.stdDev print 'Histogram:', impstats.histogram  Note that the field 'histogram' returns a list. Note also that fields labeled "protected" cannot be accessed. We now know the minimum and the mean pixel values, so we select the region bounded by these lower and upper values. We do intensity threshold. ip.setThreshold(ipmin, ipmean, ImageProcessor.RED_LUT)  setThreshold is a method of class ImageProcessor and has three arguments. 1. lower threshold value 2. upper threshold value 3. LUT update (a number) The first and the second arguments are obvious what they are, but the third one probably is not. This value determines the highlighting color of thresholded area. There are three different ways: 1. red (RED_LUT, 0) 2. black and white (BLACK_AND_WHITE_LUT, 1) 3. no update (NO_LUT_UPDATE, 2) 4. green and blue (OVER_UNDER_LUT, 3) So if you want red highlighted area, then you could put 0 as the third argument. How could we know such options are there? They are actually listed as constant field values: http://rsbweb.nih.gov/ij/developer/api/constant-values.html#ij.process.ImageProcessor.RED_LUT They are a special type of field values associated with class ImageProcessor and hold constant values. For example, ImageProcessor.RED_LUT is 0. Check this by running the following code. print ImageProcessor.RED_LUT print ImageProcessor.BLACK_AND_WHITE_LUT print ImageProcessor.NO_LUT_UPDATE print ImageProcessor.OVER_UNDER_LUT  The key point here is that these field values have a very descriptive name of what they are doing. This allows the programmer to easily select one of the options of how highlighting will appear as the result of intensity thresholding. In our case we do not care because we don't need to check the thresholded area visually. However, only because the setThreshold method needs a third argument, we set it to a default value. We now measure the pixel values only those which were selected by this thresholding. For measurement, we now use [the getSatistics method of ImageStatistics class](http://rsbweb.nih.gov/ij/developer/api/ij/process/ImageStatistics.html#getStatistics(ij.process.ImageProcessor, int, ij.measure.Calibration)): static ImageStatistics getStatistics(ImageProcessor ip, int mOptions, Calibration cal) From this, we understand that • This is a static method so we can use it directly without constructing an instance. • It returns an ImageSatistics object. • There are three arguments: 1. An image processor object 2. measurement options, a number 3. Calibration object The first one is obvious, the image that we are working on. The second is a specific number that sets the measurement option. Looking through the Javadoc of ImageStatistics class, you will find a small table of "Fields inherited from interface ij.measure.Measurements". They are listed in the constant field value table. http://rsbweb.nih.gov/ij/developer/api/constant-values.html#ij.measure.Measurements.AREA As we want to get a mean value of thresholded area, we need to have the option ImageStatistics.MEAN At the same time, we also need to set the statistics to be computed only from thresholded area so ImageStatistics.LIMIT We add these values to get a specific number that set these options active. #set measurement option measOpt = ImageStatistics.MEAN + ImageStatistics.LIMIT #get statistics of thresholded area impstats = ImageStatistics.getStatistics(ip, measOpt, None)  As we do not need calibration, we set it to null (in Java) or None (in Python). Now the measurement is done so we can use the measured mean intensity to subtract it from image (background subtraction, finally) backlevel = impstats.mean ip.resetThreshold() # subtract background level ip.subtract(backlevel) print imp.getTitle(), " : background intensity - ", backlevel  Taking all these fragments into one, we now have a script that does background subtraction. import os #path information root = '/g/data/bio-it_centres_course/data/course' filepm = '--W00005--P00001--Z00000--T00000--pm-647.tif' filevsvg = '--W00005--P00001--Z00000--T00000--vsvg-sfp.tif' fullpath = os.path.join(root, filepm) #Load an image imp = IJ.openImage(fullpath) imp.show() #Get ImageProcessor object of current image ip = imp.getProcessor() #Get ImageStatistics object of current ImageProcessor Object impstats = ip.getStatistics() #Get rtatistics results ipmin = impstats.min ipmean = impstats.mean #lower threshold = minimal intensity #upper threshold = mean intensity ip.setThreshold(ipmin, ipmean, ImageProcessor.RED_LUT) #set measurement option measOpt = ImageStatistics.MEAN + ImageStatistics.LIMIT #get statistics of thresholded area impstats = ImageStatistics.getStatistics(ip, measOpt, None) backlevel = impstats.mean ip.resetThreshold() # subtract background level ip.subtract(backlevel) print imp.getTitle(), " : background intensity - ", backlevel  In this code, only plasma membrane image is processed. To do the background subtraction for VSVG image, we could duplicate the code and replace the path for "IJ.openImage(path)", but a smarter way of this is to create a function like below. import os root = '/Volumes/data/bio-it_centres_course/data/course' filepm = '--W00005--P00001--Z00000--T00000--pm-647.tif' filevsvg = '--W00005--P00001--Z00000--T00000--vsvg-cfp.tif' def backSub(imp): imp.show() ip = imp.getProcessor() impstats = ip.getStatistics() ipmin = impstats.min ipmean = impstats.mean #ipmax = impstats.max ip.setThreshold(ipmin, ipmean, ImageProcessor.RED_LUT) measOpt = ImageStatistics.MEAN + ImageStatistics.LIMIT impstats = ImageStatistics.getStatistics(ip, measOpt, None) backlevel = impstats.mean ip.resetThreshold() ip.subtract(backlevel) print imp.getTitle(), " : background intensity - ", backlevel fullpath = os.path.join(root, filepm) imp = IJ.openImage(fullpath) backSub(imp) fullpath = os.path.join(root, filevsvg) imp = IJ.openImage(fullpath) backSub(imp)  ### Nucleus segmentation A script for processing single image dataset is: https://gist.github.com/miura/9765308 from ij.plugin.filter import GaussianBlur from fiji.threshold import Auto_Local_Threshold as ALT from ij.plugin.filter import ParticleAnalyzer as PA from org.jfree.data.statistics import BoxAndWhiskerCalculator from java.util import ArrayList, Arrays import os # size of juxtanuclear region. In pixels. RIMSIZE = 15 # image background is expected to be black. Prefs.blackBackground = True # verbose output VERBOSE = False root = '/Volumes/data/bio-it_centres_course/data/course' filedapi = '--W00005--P00001--Z00000--T00000--dapi.tif' nucpath = os.path.join(root, filedapi) impN = IJ.openImage(nucpath) class InstBWC(BoxAndWhiskerCalculator): def __init__(self): pass def getOutlierBound(rt): """ Analyzes the results of the 1st partcile analysis. Since the dilation of nuclear perimeter often causes overlap of neighboring neculeus 'terrirories', such nucleus are discarded from the measurements. Small nucelei are already removed, but since rejection of nuclei depends on standard outlier detection method, outliers in both smaller and larger sizes are discarded. """ area = rt.getColumn(rt.getColumnIndex('Area')) circ = rt.getColumn(rt.getColumnIndex("Circ.")) arealist = ArrayList(Arrays.asList(area.tolist())) circlist = ArrayList(Arrays.asList(circ.tolist())) bwc = InstBWC() ans = bwc.calculateBoxAndWhiskerStatistics(arealist) #anscirc = bwc.calculateBoxAndWhiskerStatistics(circlist) if (VERBOSE): print ans.toString() print ans.getOutliers() q1 = ans.getQ1() q3 = ans.getQ3() intrange = q3 - q1 outlier_offset = intrange * 1.5 return q1, q3, outlier_offset def nucleusSegmentation(imp2): """ Segmentation of nucleus image. Nucleus are selected that: 1. No overlapping with dilated regions 2. close to circular shape. Deformed nuclei are rejected. Outputs a binary image. """ #Convert to 8bit ImageConverter(imp2).convertToGray8() #blur slightly using Gaussian Blur radius = 2.0 accuracy = 0.01 GaussianBlur().blurGaussian( imp2.getProcessor(), radius, radius, accuracy) # Auto Local Thresholding imps = ALT().exec(imp2, "Bernsen", 15, 0, 0, True) imp2 = imps[0] #ParticleAnalysis 0: prefiltering by size and circularity rt = ResultsTable() paOpt = PA.CLEAR_WORKSHEET +\ PA.SHOW_MASKS +\ PA.EXCLUDE_EDGE_PARTICLES +\ PA.INCLUDE_HOLES #+ \ # PA.SHOW_RESULTS measOpt = PA.AREA + PA.STD_DEV + PA.SHAPE_DESCRIPTORS + PA.INTEGRATED_DENSITY MINSIZE = 20 MAXSIZE = 10000 pa0 = PA(paOpt, measOpt, rt, MINSIZE, MAXSIZE, 0.8, 1.0) pa0.setHideOutputImage(True) pa0.analyze(imp2) imp2 = pa0.getOutputImage() # Overwrite imp2.getProcessor().invertLut() impNuc = imp2.duplicate() ## for the ring. #impNuc = Duplicator().run(imp2) #Dilate the Nucleus Area ## this should be 40 pixels in Cihan's method, but should be smaller. rf = RankFilters() rf.rank(imp2.getProcessor(), RIMSIZE, RankFilters.MAX) #Particle Analysis 1: get distribution of sizes. paOpt = PA.CLEAR_WORKSHEET +\ PA.SHOW_NONE +\ PA.EXCLUDE_EDGE_PARTICLES +\ PA.INCLUDE_HOLES #+ \ # PA.SHOW_RESULTS measOpt = PA.AREA + PA.STD_DEV + PA.SHAPE_DESCRIPTORS + PA.INTEGRATED_DENSITY rt1 = ResultsTable() MINSIZE = 20 MAXSIZE = 10000 pa = PA(paOpt, measOpt, rt1, MINSIZE, MAXSIZE) pa.analyze(imp2) #rt.show('after PA 1') #particle Analysis 2: filter nucleus by size and circularity. #print rt1.getHeadings() if (rt1.getColumnIndex('Area') > -1): q1, q3, outlier_offset = getOutlierBound(rt1) else: q1 = MINSIZE q3 = MAXSIZE outlier_offset = 0 print imp2.getTitle(), ": no Nucleus segmented,probably too many overlaps" paOpt = PA.CLEAR_WORKSHEET +\ PA.SHOW_MASKS +\ PA.EXCLUDE_EDGE_PARTICLES +\ PA.INCLUDE_HOLES #+ \ # PA.SHOW_RESULTS rt2 = ResultsTable() pa = PA(paOpt, measOpt, rt2, q1-outlier_offset, q3+outlier_offset, 0.8, 1.0) pa.setHideOutputImage(True) pa.analyze(imp2) impDilatedNuc = pa.getOutputImage() #filter original nucleus filteredNuc = ImageCalculator().run("AND create", impDilatedNuc, impNuc) return filteredNuc imp2 = impN.duplicate() impfilteredNuc = nucleusSegmentation(imp2) impfilteredNuc.show()  https://gist.github.com/miura/9644827 ### ROI generator, Measurements & Transport Ratio • Accessing ResultsTable, read and write • ParticleAnalysis class • from x import y as z • ROI, ShapeRoi classes • Implementing Java Interface • ImageCalculator class import os, sys from ij.plugin.filter import ParticleAnalyzer as PA # size of juxtanuclear region. In pixels. RIMSIZE = 15 def roiRingGenerator(r1): """ Create a band of ROI outside the argument ROI. See Liebel (2003) Fig. 1 """ #r1 = imp.getRoi() r2 = RoiEnlarger.enlarge(r1, RIMSIZE) sr1 = ShapeRoi(r1) sr2 = ShapeRoi(r2) return sr2.not(sr1) def roiEnlarger(r1): """ Enlarges ROI by a defined iterations. """ return ShapeRoi(RoiEnlarger.enlarge(r1, RIMSIZE)) def measureROIs(imp, measOpt, thisrt, roiA, backint, doGLCM): """ Cell-wise measurment using ROI array. """ analObj = Analyzer(imp, measOpt, thisrt) for index, r in enumerate(roiA): imp.deleteRoi() imp.setRoi(r) analObj.measure() maxint = thisrt.getValue('Max', thisrt.getCounter()-1) saturation = 0 if ( maxint + backint) >= 4095: saturation = 1 if (VERBOSE): print 'cell index ', index, 'maxint=', maxint thisrt.setValue('CellIndex', thisrt.getCounter()-1, index) thisrt.setValue('Saturation', thisrt.getCounter()-1, saturation) if (doGLCM): imp.deleteRoi() #measureTexture(imp, thisrt, roiA) ## loading files root = '/Volumes/data/bio-it_centres_course/data/course' fileNuc = 'dapiSegmented.tif' filepm = 'pm-647BackSub.tif' filevsvg = 'vsvg-cfpBackSub.tif' nucfull = os.path.join(root, fileNuc) pmfull = os.path.join(root, filepm) vsvgfull = os.path.join(root, filevsvg) impfilteredNuc = IJ.openImage(nucfull) impPM = IJ.openImage(pmfull) impVSVG = IJ.openImage(vsvgfull) wnumber = '00005' rtallcellPM = ResultsTable() rtjnucVSVG = ResultsTable() rtallcellVSVG = ResultsTable() intmax = impfilteredNuc.getProcessor().getMax() if intmax == 0: #return rtallcellPM, rtjnucVSVG, rtallcellVSVG exit() impfilteredNuc.getProcessor().setThreshold(1, intmax, ImageProcessor.NO_LUT_UPDATE) nucroi = ThresholdToSelection().convert(impfilteredNuc.getProcessor()) nucroiA = ShapeRoi(nucroi).getRois() #print nucroiA allcellA = [roiEnlarger(r) for r in nucroiA] jnucroiA = [roiRingGenerator(r) for r in nucroiA] measOpt = PA.AREA + PA.MEAN + PA.CENTROID + PA.STD_DEV + PA.SHAPE_DESCRIPTORS + PA.INTEGRATED_DENSITY + PA.MIN_MAX +\ PA.SKEWNESS + PA.KURTOSIS + PA.MEDIAN + PA.MODE ## All Cell Plasma Membrane intensity measureROIs(impPM, measOpt, rtallcellPM, allcellA, 0, True) meanInt_Cell = rtallcellPM.getColumn(rtallcellPM.getColumnIndex('Mean')) print "Results Table rownumber:", len(meanInt_Cell) # JuxtaNuclear VSVG intensity measureROIs(impVSVG, measOpt, rtjnucVSVG, jnucroiA, 0, False) meanInt_jnuc = rtjnucVSVG.getColumn(rtjnucVSVG.getColumnIndex('Mean')) # AllCell VSVG intensity measureROIs(impVSVG, measOpt, rtallcellVSVG, allcellA, 0, True) meanInt_vsvgall = rtallcellVSVG.getColumn(rtallcellVSVG.getColumnIndex('Mean')) for i in range(len(meanInt_Cell)): if meanInt_Cell[i] != 0.0: transportR = meanInt_jnuc[i] / meanInt_Cell[i] transportRall = meanInt_vsvgall[i] / meanInt_Cell[i] else: transportR = float('inf') transportRall = float('inf') rtjnucVSVG.setValue('TransportRatio', i, transportR) rtallcellVSVG.setValue('TransportRatio', i, transportRall) rtjnucVSVG.setValue('WellNumber', i, int(wnumber)) rtallcellVSVG.setValue('WellNumber', i, int(wnumber)) rtallcellPM.setValue('WellNumber', i, int(wnumber)) rtallcellVSVG.show('AllCell')  ### Pre-screening of Images. We will not go over the details of full processing, but the pre-screening code is a good example of accessing multiple files and process images. https://gist.github.com/miura/9684453 Pre-screen results are saved under /g/data/bio-it_centres_course/data/VSVG/prescreen  ## Final Code Assembling all the fragments as a single script, here is the final code to submit. https://github.com/miura/HTManalysisCourse/blob/master/measTransportBatch3.py An important command that appears only here is 'sys.argv, which holds the argument added in the command line. Try the following to see what is happening. Create a file 'testargv.py' with following code. import sys print 'The argument given: ', sys.argv  Save the file, then run it from the command line. fiji --mem1000m testargv.py testarg1 testarg2  The output should look like miura@compute-n026<19> [/g/data/bio-it_centres_course/data/course ] fiji --mem=2000m testargv.py testarg1 testarg2 No GUI detected. Falling back to headless mode. No GUI detected. Falling back to headless mode. the argument givien: ['testargv.py', 'testarg1', 'testarg2']  Within the running Jython instance, command line arguments are passed to the program via sys.argv. In our case, we pass name of the plate to process to the script. As an example, if you want to process the first plate (0001-03--2005-08-01) and save outputs (measurements) to a folder named 'out_kota', the command will be fiji --mem2000m 0001-03--2005-08-01 out_kota  As there is no line in the script to create a new folder, 'out_kota' should be created manually before doing this. If there is any problem with memory, try increasing it a bit, such as --mem2500m or so. ## Generating Job Array Here is the code written by Clemens. https://github.com/miura/HTManalysisCourse/blob/master/fiji_script.sh #!/bin/bash -e # Read-only filesystem error occurs sometimes # Appears to be due to mounting issues or file system limitations # Retry if it occurs BASE_DIR="/g/data/bio-it_centres_course/data/VSVG" SUB_DIR="${BASE_DIR}/fiji-sub"
LOG_DIR="${BASE_DIR}/fiji-log" PRE_DIR="${BASE_DIR}/prescreen"
# The first version of the python script gave an error for some plates
#JYTHON_SCRIPT="/g/almf/software/scripts2/measTransportBatch.py"
#JYTHON_SCRIPT="/g/almf/software/scripts2/measTransportBatch2.py"
JYTHON_SCRIPT="/g/almf/software/scripts2/measTransportBatch3.py"

cd ${BASE_DIR} [ -d fiji-sub ] || mkdir fiji-sub [ -d fiji-log ] || mkdir fiji-log [ -d fiji-out ] || mkdir fiji-out for a in ls${BASE_DIR} | grep "\-\-"
do
PLATE=${a} # Check if prescreen exists if [ ! -f${PRE_DIR}/${PLATE}.csv ]; then echo "Skipping plate${PLATE}: no pre-screen data."
continue
fi
OO="${LOG_DIR}/${PLATE}-out.txt"
EO="${LOG_DIR}/${PLATE}-err.txt"
TARGET_DIR="${PLATE}" OUTPUT_DIR="fiji-out" # Continue if output already exists - no need to re-run if [ -f${OUTPUT_DIR}/${PLATE}--PMall.csv ]; then echo "Skipping plate${PLATE}: results already exist."
continue
fi

# Write the bsub script
echo "#!/bin/bash
#BSUB -oo \"${OO}\" #BSUB -eo \"${EO}\"
#BSUB -M 8000
#BSUB -R select[mem>8000] -R rusage[mem=8000]
ulimit -c 0
export PATH=/g/almf/software/bin2:\${PATH} echo \"The job started:\" java -version JYTHON_SCRIPT=${JYTHON_SCRIPT}
TARGET_DIR=${TARGET_DIR} OUTPUT_DIR=${OUTPUT_DIR}
echo \"Fiji started:\"
fiji --headless --mem=1000m ${JYTHON_SCRIPT}${TARGET_DIR} ${OUTPUT_DIR}" >${SUB_DIR}/${PLATE}.bsub # Submit the job to the cluster chmod +x${SUB_DIR}/${PLATE}.bsub bsub <${SUB_DIR}/\${PLATE}.bsub
done

`
documents/140326jointcourse.txt · Last modified: 2017/06/07 03:08 by kota