A minimal use of opencv: toward chromosome shape analysis.

In an attempt to analyse the chromosomes shape, I need to access to the coordinates of points belonging to the contour or to the skeleton of chromosomes. Opencv is a rich vision library, usable from python. Here, I wrote a minimal script to access the pixels of the chromosome contour:

part15.png:image used to test opencv

import cv
definition of some colors
_red =  (0, 0, 255, 0);
_green =  (0, 255, 0, 0);
_white = cv.RealScalar (255)
_black = cv.RealScalar (0)
# create window and display the original picture in it
path="/home/claire/Applications/ProjetPython/test opencv/part15.png"
image = cv.LoadImage(path,cv.CV_LOAD_IMAGE_GRAYSCALE)
binary=image>0
cv.NamedWindow ("image", 1)
cv.ShowImage ("image", image)
# create the storage area
storage = cv.CreateMemStorage (0)
# find the contours
contours = cv.FindContours(image,
                               storage,
                               cv.CV_RETR_TREE,
                               cv.CV_CHAIN_APPROX_SIMPLE,
                               (0,0))
cv.DrawContours(image,contours,_red,_black,max_level=2)
print len(contours)
print contours[0],contours[1],contours[61]
cv.WaitKey (0)

• The first thing I noticed is that without cv.WaitKey(0), the image is not visible, it is displayed and erased immediately.
• The object contours is a list of pixels coordinates, I didn't check yet if it is an ordered list of pixels.

My feeling about opencv is that it is a world itself, may be not mixing very well with other library.

Making a MFISH image with scikits.image

I got a great help from the scikits.image google group to make a MFISH image but with a nicer code than the previous one.

Making a MFISH image

The first trial with numpy was not very successful. After code simplification, I get a result:

RGB image combining five monospectral images

The code is:

# -*- coding: utf-8 -*-
"""
Created on Mon May 23 14:07:40 2011

@author: jean-pat
"""

import numpy as np
import os
import pylab
#from scipy import linalg as la
#from scipy import ndimage as nd
import scipy as sp
user=os.path.expanduser("~")
workdir=os.path.join(user,"Applications","ImagesTest","MFISH")
blue="Aqua.tif"
green="Green.tif"
gold="Gold.tif"
red="Red.tif"
frd="FarRed.tif"
complete_path=os.path.join(workdir,blue)
i1=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,green)
i2=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,gold)
i3=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,red)
i4=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,frd)
i5=sp.misc.imread(complete_path)

R=0.8*i5+0.15*i4+0.05*i3+0.00*i2+0.00*i1
V=0.00*i5+0.05*i4+0.05*i3+0.80*i2+0.10*i1
B=0.00*i5+0.00*i4+0.05*i3+0.15*i2+0.80*i1
shape=i5.shape
rgb = np.zeros((shape[0],shape[1],3),dtype=np.uint8)
rgb[:,:,0]=np.copy(R)
rgb[:,:,1]=np.copy(V)
rgb[:,:,2]=np.copy(B)
pylab.subplot(111,frameon=False, xticks=[], yticks=[])
pylab.imshow(rgb)
pylab.show()

I received new advices today

# -*- coding: utf-8 -*-
"""
Created on Mon May 23 14:07:40 2011

@author: jean-pat
"""

import numpy as np
import os
import pylab
#from scipy import linalg as la
#from scipy import ndimage as nd
import scipy as sp
user=os.path.expanduser("~")
workdir=os.path.join(user,"Applications","ImagesTest","MFISH")
blue="Aqua.tif"
green="Green.tif"
gold="Gold.tif"
red="Red.tif"
frd="FarRed.tif"
complete_path=os.path.join(workdir,blue)
i1=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,green)
i2=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,gold)
i3=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,red)
i4=sp.misc.imread(complete_path)
#
complete_path=os.path.join(workdir,frd)
i5=sp.misc.imread(complete_path)
######################
multi_spectr = np.dstack([i1,i2,i3,i4,i5])
conv_matrix = np.array([
[0.8,0.15,0.05,0,0],
[0,0.05,0.05,0.8,0.1],
[0,0,0.05,0.15,0.8]
])
 
shape = multi_spectr.shape
multi_spectr = multi_spectr.reshape(shape[0]*shape[1],shape[2])
rgb = np.dot(multi_spectr,conv_matrix.T).reshape(shape[:2]+(3,))
rgb8=np.uint8(rgb)
########################
#R=0.8*i5+0.15*i4+0.05*i3+0.00*i2+0.00*i1
#V=0.00*i5+0.05*i4+0.05*i3+0.80*i2+0.10*i1
#B=0.00*i5+0.00*i4+0.05*i3+0.15*i2+0.80*i1
#shape=i5.shape
#rgb = np.zeros((shape[0],shape[1],3),dtype=np.uint8)
#rgb[:,:,0]=np.copy(R)
#rgb[:,:,1]=np.copy(V)
#rgb[:,:,2]=np.copy(B)
pylab.subplot(111,frameon=False, xticks=[], yticks=[])
pylab.imshow(rgb8)
pylab.show()

That yields:

MFISH image to RGB with image5D

Wade Schwartzkopf published a MFISH dataset  in 2000. Image5D (a plugin for imagej) can display five different fluorochromes simultaneously:

Left:raw png images. Top right: RGB image done with Image5D. Bottom right: "karyotype" image with false colors. 

Far red

Red

Spectrum gold

Spectrum green

Spectrum aqua

Dapi

After chromosome classification, a karyotype image can be produced. This is a grayscaled image associated to a  LUT, to be displayed with false colors:

Image5D can export a rgb image from the five channels image:

RGB image generated after background substraction of the five raw images

Highlighting Interstitial Telomeric Spots in interphasic cells (HITS-FISH)

Telomeric fusions underscore telomeric binding proteins disfunctions. Fusions occur in senescent cells between short telomeres, they are detected in metaphasic cells by dicentric chromosomes most generally with no telomeric spots at the fusion point.

In some SV40 transformed fibroblast cell line, interstitial telomeric spots (ITS) are as bright as terminal telomeric spots, this indicates some uncoupling between telomere lenght and their capacity to block telomeric fusions. 

ITS observed in a precrisis human SV40 transformed cell line (Ducray et al., oncogene 1999)

From metaphasic chromosomes observation, fusions retaining visible telomeric spots are rare events but may be more frequent in interphasic cells.

ITS from a SV40 transformed cell line

Bal31 in-situ:

Several nucleases or DNA polymerases  were used in assays requiring fixed material (chromosome banding+restriction endonuclease, in-situ nick translation, fluorescent labelling with TUNEL, in-situ PCR). The exonuclease Bal31 removes totally terminal telomeric sequences of naked DNA after some hours of digestion (R C Allshire, M Dempster, and N D Hastie, NAR 1989; de Lange et al, MCB 1990) . In-situ,  the terminal telomeric sequences, should be removed by Bal31 digestion.  If telomeric spots can be detected in nuclei after Bal31 digestion, those spots should be  ITS. Monitoring ITS may be a tool of interest  to probe telesome functions at the scale of unique cell, or to monitor the state of a tissue.

Pig cells as model cells?

On human metaphasic chromosomes, a bright ITS indicates an  abnormal cell but for cheaper DAPI stained chromosomes rearrangements too. Normal pig karyotype  shows interstitial telomeric spots:

Normal diploid pig metaphase with a chromosomes pair having ITS.

After QFISH on interphasic pig nuclei, there are two ITS of same magnitude than terminal spots, these two pig ITS remain hidden in the forest of terminal telomeric spots:

  

Pig nuclei (25x magnification)

Some conditions of Bal31 digestion should retain these two ITS and remove totally terminal telomeric spots yielding the possibility to observe telomeric fusions in nuclei.

Assay:

Tuning the conditions of Bal31 digestion is a prerequisite step. Supposing that the exonuclease Bal31 can remove terminal sequence from fixed material, when this enzyme is usually used to digest purified DNA in solution. 

If pig cells are not available, performing Bal31 digestion on control human cells such lymphocytes should tell if Bal31 can remove telomeres in-situ.

From conventionnal cytogenetic preparation of pig cells (fibroblasts,lymphocytes), digestion time, numbers of Bal31 units could be set up as follow : 

Possible protocol for HITS FISH

The protocol modified from Lansdorp et al. (Hum Mol Genet. 1996 May;5(5):685-91.) could be also adapted to liquid hybridization for flow FISH. Under optimal conditions pig cells could be used as internal control for Bal 31 digestion by mixing human and pig cells in the same cytogenetic preparation (9:1 ratio).

Human cells can be then distinguished from pig cells with an additionnal oligo probe. In the following example (25x magnification), simultaneous hybridization with Cy3-PNA-(CCCTAA)3 (displayed in green) and Cy5-DNA (displayed in red) oligo probe was performed (unpublished data):

Candidate cells to test the assay, may be SV40 precrisis cells with highly rearranged chromosomes and numerous dicentrics, or better cells with controlled mutations.

A better chromosomes gallery with size filtering

In a previous post, I tried to code using pylab/matplotlib an image gallery to get a result similar to make montage in imageJ.
I have written a python function called makeMosaic. It take three arguments:

• a list of particules  consisting in monospectral images (a numpy array)
• a minimal size (pixel)
• a maximal size

In a python script the function is used like that:

makeMosaic(Particles,100,2000)

Two examples of result implemented in a script:

Gallery in the range 100-10000

Gallery in the range 100-2000

Original DAPI stained metaphase

Particle : a python class to compute particle anisotropy by morphological erosion

A python class called particle, is implemented to handle the chromosomes extracted from a metaphase image:

Original segmented chromosomes

The class particle  is instanciated with the counterstain image (DAPI), it is possible to :

1. display a particle vertically
2. compute the convex hull area. The convex hull area is computed with from the coordinates of the convex hull vertices.
3. produce the morphological rose des directions, which describe the particle anisotropy.

Blue curve:orientation curve (rose des direction) for the original particle image (top left). Green curve: orientation curve after high pass filtering (top right).

The script written to test the class, instanciates two particle objects, for example:

•  p1=particle(im), where im is a DAPI image of segmented chromosomes.

It is then possible to get the morphological "rose des directions", that is the number of morphological erosion by a 1x3 structuring element, necessary to erode totally the particle in a given direction:

• theta1,rosedirections,VImage=p1.orientationByErosion(10) 

• The parameter set to 10 is the rotation step
• theta: the particle principal orientation
• rosedirections: a numpy array to plot the erosion number vs the rotation angle
• VImage: rotated particle image

 The particle convex hull is obtained with:

• p1.cvxhull_area()

The script must be modified to fit  the image path.

Numpy, scipy, matplotlib, pymorph, readmagick, mahotas must be installed. However since original implementation. Mahotas itself or scikit image could be used instead of readmagick to load the image.

With two features, as the particle area and the ratio between the particle area and the convex hull area, it should be possible to start to classify the particles into four categories:

• non overlapping chromosomes,
• overlapping chromosomes,
• nuclei, 
• remaining small stuffs .

Code on Gist 

anisotropy_RoseDirections.py

Ipython notebook here

http://nbviewer.ipython.org/gist/jeanpat/7673858