When looking at the microscopy images of ER exit sites we noticed that after a certain perturbation they change the distribution around the nuclei. ER exit sites in normal cells are more concentrated on one side. After the perturbation, the ER exit sites uniformly distribute around the nuclei. We wanted to quantify this phenotype, so I developed a simple algorithm to express the uniformity of ER exit sites distribution with a single number.
A theoretic example
Generate three different samples of angles with mean angle 0. I use von Mises distribution, which is a continuous probability distribution on the unit circle. It may be considered as a circular analogue of the normal distribution. Then calculate an 8-bin, normalized histogram of angles. The Polar Distribution Score is calculated as the difference the sample of angles and uniform distribution. The difference is expressed as the absolute sum of histogram differences.
import numpy as np
import matplotlib.pyplot as plt
fig,ax = plt.subplots(2,3,figsize=(15,10),facecolor="w")
# define histogram bins and radius
bins=np.linspace(-np.pi,np.pi,9)
r=np.random.uniform(0.2,0.8,100)
# Iterate over populations
for i,k in enumerate([0,1,10]):
# Angles are sampled from von Mises distribution:
# mu - mean angle, kappa - spread (>=0), size - sample size
theta=np.random.vonmises(mu=0, kappa=k, size=100)
# Plot on polar projection.
ax[0,i].remove()
ax[0,i] = fig.add_subplot(2, 3,i+1, projection = 'polar')
c = ax[0,i].scatter(theta, r, alpha=0.75)
ax[0,i].set_ylim((0,1))
ax[0,i].set_title(f"Mean angle 0, spread {k}")
# Plot histogram
h = ax[1,i].hist(theta, bins=bins, density=True)
ax[1,i].set_ylim((0,1))
# Bin density of uniform distribution is 1/(8 * bin_width)
bin_width = np.diff(bins)[0]
uniform_hist = 1/(8 * bin_width)
ax[1,i].axhline(uniform_hist, c="r")
# PDS is the absolute sum of differences
pds = np.sum(np.abs(h[0]-uniform_hist))
ax[1,i].text(0, 0.9, f"PDS: {pds:.2f}", ha="center")
A practical demonstration
Let’s demonstrate this on a real life example. I’ll take images above, segment the cells and ER exit sites. Then I’ll convert the cartesian coordinates to polar and calculate the PDS.
from skimage import io
from skimage.filters import threshold_local, gaussian, threshold_otsu
from skimage.segmentation import clear_border
from skimage.measure import label, regionprops_table
import pandas as pd
img_disrupted_eres=io.imread("pds-cells-rhob-green.png")
img_normal_eres=io.imread("pds-cells-neg-green.png")
img_disrupted_eres_blurred=gaussian(img_disrupted_eres,80)
img_normal_eres_blurred=gaussian(img_normal_eres,80)
cells_disrupted_mask=img_disrupted_eres_blurred>threshold_local(img_disrupted_eres_blurred,501,offset=0.001)
cells_normal_mask=img_normal_eres_blurred>threshold_local(img_normal_eres_blurred,501,offset=0.001)
#Eres
eres_disrupted_mask=img_disrupted_eres>threshold_otsu(img_disrupted_eres)
eres_normal_mask=img_normal_eres>threshold_otsu(img_normal_eres)
eres_disrupted=label(eres_disrupted_mask)
eres_normal=label(eres_normal_mask)
#Nuclei
img_disrupted_nuc=io.imread("pds-cells-rhob-blue.png")
img_normal_nuc=io.imread("pds-cells-neg-blue.png")
img_disrupted_nuc_blurred=gaussian(img_disrupted_nuc,15)
img_normal_nuc_blurred=gaussian(img_normal_nuc,15)
nuc_disrupted_mask=clear_border(img_disrupted_nuc_blurred>threshold_local(img_disrupted_nuc_blurred,501,offset=0.001))
nuc_normal_mask=clear_border(img_normal_nuc_blurred>threshold_local(img_normal_nuc_blurred,501,offset=0.001))
# Link cells and nuclei
cells_normal_mask=clear_border(np.logical_or(nuc_normal_mask,cells_normal_mask))
cells_disrupted_mask=clear_border(np.logical_or(nuc_disrupted_mask,cells_disrupted_mask))
cells_disrupted=label(cells_disrupted_mask)
cells_normal=label(cells_normal_mask)
nuc_disrupted=nuc_disrupted_mask*cells_disrupted
nuc_normal=nuc_normal_mask*cells_normal
# Get region properties for nuclei.
nuc_disrupted_rp=pd.DataFrame(regionprops_table(nuc_disrupted,properties=["label","centroid"])).set_index("label")
nuc_normal_rp=pd.DataFrame(regionprops_table(nuc_normal,properties=["label","centroid"])).set_index("label")
# Link ERES to cells and get region props.
eres_disrupted_rp=[]
for c in np.unique(cells_disrupted[cells_disrupted>0]):
mask=eres_disrupted*(cells_disrupted==c)
eres_rp=pd.DataFrame(regionprops_table(mask,
properties=["label","centroid"]))
eres_rp["cell"]=c
eres_disrupted_rp.append(eres_rp)
eres_disrupted_rp=pd.concat(eres_disrupted_rp).set_index("label")
eres_normal_rp=[]
for c in np.unique(cells_normal[cells_normal>0]):
mask=eres_normal*(cells_normal==c)
eres_rp=pd.DataFrame(regionprops_table(mask,
properties=["label","centroid"]))
eres_rp["cell"]=c
eres_normal_rp.append(eres_rp)
eres_normal_rp=pd.concat(eres_normal_rp).set_index("label")
# Subtract nuceli centroids to center the coordinates
eres_disrupted_rp[["centroid_centered-0","centroid_centered1"]]=eres_disrupted_rp.apply(lambda x: x[["centroid-0","centroid-1"]]-
nuc_disrupted_rp.loc[x["cell"],["centroid-0","centroid-1"]],axis=1)
eres_normal_rp[["centroid_centered-0","centroid_centered1"]]=eres_normal_rp.apply(lambda x: x[["centroid-0","centroid-1"]]-
nuc_normal_rp.loc[x["cell"],["centroid-0","centroid-1"]],axis=1)
# Convert to polar coordinates
eres_disrupted_rp["rho"] =
np.sqrt(eres_disrupted_rp["centroid_centered-1"]**2 +
eres_disrupted_rp["centroid_centered-0"]**2)
eres_disrupted_rp["phi"] =
np.arctan2(eres_disrupted_rp["centroid_centered-0"],
eres_disrupted_rp["centroid_centered-1"])
eres_normal_rp["rho"] = np.sqrt(eres_normal_rp["centroid_centered1"]**2 + eres_normal_rp["centroid_centered-0"]**2)
eres_normal_rp["phi"] = np.arctan2(eres_normal_rp["centroid_centered0"], eres_normal_rp["centroid_centered-1"])
# zero mean angle normalization
eres_disrupted_rp["phi_norm"]=(eres_disrupted_rp.groupby("cell").apply (lambda x: x["phi"]-np.angle(np.sum(np.exp([complex(0,f) for f in x["phi"]])))).to_numpy()+2*np.pi)%(2*np.pi)
eres_normal_rp["phi_norm"]=(eres_normal_rp.groupby("cell").apply(lambda x: x["phi"]-np.angle(np.sum(np.exp([complex(0,f) for f in x["phi"]])))).to_numpy()+2*np.pi)%(2*np.pi)
# histogram and pds per cell
def calculate_pds(a):
bins=np.linspace(0,2*np.pi,9)
bin_width = np.diff(bins)[0]
uniform_hist = 1/(8 * bin_width)
a=(a+2*np.pi)%(2*np.pi)
h=np.histogram(a,bins=bins, density=True)
pds = np.sum(np.abs(h[0]-uniform_hist))
return(pds)
eres_disrupted_pds=eres_disrupted_rp.groupby("cell").agg(pds=("phi_norm",calculate_pds))
eres_normal_pds=eres_normal_rp.groupby("cell").agg(pds=("phi_norm",calculate_pds))
# Plot
bins=np.linspace(0,2*np.pi,9)
bin_width = np.diff(bins)[0]
uniform_hist = 1/(8 * bin_width)
ncells=np.maximum(len(eres_disrupted_rp["cell"].unique()),len(eres_normal_rp["cell"].unique()))
fig,ax = plt.subplots(4,ncells,figsize=(5*ncells,22),facecolor="w")
# Iterate over disrupted cells
for i,k in enumerate(eres_disrupted_rp["cell"].unique()):
phi=eres_disrupted_rp.loc[eres_disrupted_rp["cell"]==k,"phi_norm"].to_numpy()
rho=eres_disrupted_rp.loc[eres_disrupted_rp["cell"]==k,"rho"].to_numpy()
# Plot on polar projection.
ax[0,i].remove()
ax[0,i] = fig.add_subplot(4, ncells,i+1, projection = 'polar')
c = ax[0,i].scatter(phi, rho, alpha=0.75,s=10)
ax[0,i].set_title(f"Disrupted cell {k}")
# Plot histogram
h = ax[1,i].hist(phi, bins=bins, density=True)
ax[1,i].set_ylim((0,1))
ax[1,i].axhline(uniform_hist, c="r")
# PDS is the absolute sum of differences
ax[1,i].text(np.pi, 0.9, f"PDS:
{eres_disrupted_pds.loc[k,'pds']:.2f}", ha="center")
# Iterate over normal cells
for i,k in enumerate(eres_normal_rp["cell"].unique()):
phi=eres_normal_rp.loc[eres_normal_rp["cell"]==k,"phi_norm"].to_numpy()
rho=eres_normal_rp.loc[eres_normal_rp["cell"]==k,"rho"].to_numpy()
# Plot on polar projection.
ax[2,i].remove()
ax[2,i] = fig.add_subplot(4, ncells,i+7, projection = 'polar')
c = ax[2,i].scatter(phi, rho, alpha=0.75,s=10)
ax[2,i].set_title(f"Normal cell {k}")
# Plot histogram
h = ax[3,i].hist(phi, bins=bins, density=True)
ax[3,i].set_ylim((0,1))
ax[3,i].axhline(uniform_hist, c="r")
ax[3,i].text(np.pi, 0.9, f"PDS:{eres_normal_pds.loc[k,'pds']:.2f}", ha="center")
Disrupted cells have a lower PDS, because the distribution of ER exit sites around the nuclei is more uniform. Normal cells (rows 3 and 4 above) have higher PDS, because ER exit sites are concentrated on one edge.
The segmentation and measurements are done with scikit-image. If working with a larger dataset it would be easier to use CellProfiler to do the segmentation. Then we would import the .csv files to Jupyter, center the coordinates and do the rest.