Quantifying pixel intensities of an confocal image in Python

In this post, you are going to learn how to read an image acquired by a confocal microscope and how to measure pixel intensities. Here, I will read in an ZEISS Laser Scanning Confocal Microscope (LSM) file in Python as an example.

from pathlib import Path
import tifffile
import matplotlib.pyplot as plt
from skimage.filters import threshold_otsu

Suppose that you have a folder in which multiple image files are located. First, create a Python list containing a path to each image file. Here

data_dir = Path('path/to/folder')
paths_to_imgs = list(data_dir.glob('*.lsm'))  # Search for files end with .lsm
paths_to_imgs.sort()

Then you can read in one of the image using the path.

# Read in the first image in the folder
img = tifffile.imread(paths_to_imgs[0])

Extract axes corresponding to channel, height, and width.

img.shape
> (1, 1, 4, 2048, 2048)

img = img[0,0,:,:,:]

The image has 4 channels with 2048×2048 pixels. Display the each channel in gray-scale.

plt.figure(figsize=(16, 4))

channel_names = ['Channel 1', 'Channel 2', 'Channel 3', 'Channel 4']

for i in range(4):
    plt.subplot(1,4,i+1)
    plt.title(channel_names[i])
    plt.imshow(img[i], vmax=100, cmap='gray')  # Set vmax for enhancing contrast

Here you can Channel 2 and Channel 3 show GFP and RFP fluorescence, respectively. Now let’s compare GFP pixel intensities between the RFP-positive and RFP-negative areas.

I will create a binary mask for the RFP-positive area by thresholding.

plt.hist(img[2].flatten(), bins=256)
plt.show()

# Set a threshold
mask = img[2] >= 15

plt.figure(figsize=(12,4))
plt.subplot(131), plt.imshow(img[2], cmap='gray')
plt.subplot(132), plt.imshow(mask, cmap='gray')
plt.subplot(133), plt.imshow(~mask, cmap='gray')
plt.show()

Instead of manually choosing a threshold, you can select a threshold by the Otsu algorithm like below:

thresh = threshold_otsu(img[2])

You can calculate the mean pixel intensities of GFP signal in the RFP-positive and RFP-negative areas.

print(f'Mean pixel intensity (RFP-positive area): {img[1][mask].mean():.2f}')
print(f'Mean pixel intensity (RFP-negative area): {img[1][~mask].mean():.2f}')

The result will look like this:

Mean pixel intensity (RFP-positive area): 13.67
Mean pixel intensity (RFP-negative area): 13.65

You can write a function performing the above manipulation and apply to all the images in the folder.

def quantify(path):
    image = tifffile.imread(path)
    img = img[0,0,:,:,:]
    mask = image[2] >= 15
    GFP_RFP_neg = image[1][~mask].mean()
    GFP_RFP_pos = image[1][mask].mean()
    
    return GFP_RFP_neg, GFP_RFP_pos


# Create lists for appending the mean values of GFP intensities
GFP_intensities_RFP_neg = []
GFP_intensities_RFP_pos = []

for path in paths_to_imgs:
    GFP_RFP_neg, GFP_RFP_pos = quantify(path)
    GFP_intensities_RFP_neg.append(GFP_control)
    GFP_intensities_RFP_pos.append(GFP_experiment)