https://github.com/turleyjm/cell-division-dl-plugin
Tip revision: f6c3290670f23524d47ccddc873ef5673eda4b3c authored by turleyjm on 11 July 2024, 02:09:47 UTC
Update README.md
Update README.md
Tip revision: f6c3290
databases.py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy as sp
import skimage as sm
import skimage.io
import skimage.measure
import skimage.feature
from shapely.geometry import Polygon
import tifffile
import cellProperties as cell
import utils as util
plt.rcParams.update({"font.size": 12})
# -------------------
# This scripted re-organises and sorts data ready for analysis.
# -------------------
def dist(polygon, polygon0):
[x1, y1] = cell.centroid(polygon)
[x0, y0] = cell.centroid(polygon0)
return ((x0 - x1) ** 2 + (y0 - y1) ** 2) ** 0.5
def angleDiff(theta, phi):
diff = theta - phi
if abs(diff) > 90:
if diff > 0:
diff = 180 - diff
else:
diff = 180 + diff
return abs(diff)
def findtcj(polygon, img):
centroid = cell.centroid(polygon)
x, y = int(centroid[0]), int(centroid[1])
img = 1 - img / 255
img = np.asarray(img, "uint8")
imgLabel = sm.measure.label(img, background=0, connectivity=1)
label = imgLabel[x, y]
contour = sm.measure.find_contours(imgLabel == label, level=0)[0]
# imgLabelrc = util.imgxyrc(imgLabel)
# imgLabelrc[imgLabelrc == label] = round(1.25 * imgLabelrc.max())
# imgLabelrc = np.asarray(imgLabelrc, "uint16")
# tifffile.imwrite(f"results/imgLabel{filename}.tif", imgLabelrc)
if label == 0:
print("label == 0")
zeros = np.zeros([512, 512])
zeros[imgLabel == label] = 1
for con in contour:
zeros[int(con[0]), int(con[1])] = 1
struct2 = sp.ndimage.generate_binary_structure(2, 2)
dilation = sp.ndimage.morphology.binary_dilation(zeros, structure=struct2).astype(
zeros.dtype
)
dilation[zeros == 1] = 0
# dilationrc = util.imgxyrc(dilation)
# dilationrc = np.asarray(dilationrc, "uint16")
# tifffile.imwrite(f"results/dilation{filename}.tif", dilationrc)
tcj = np.zeros([512, 512])
diff = img - dilation
tcj[diff == -1] = 1
tcj[tcj != 1] = 0
outerTCJ = skimage.feature.peak_local_max(tcj)
tcjrc = util.imgxyrc(tcj)
tcjrc = np.asarray(tcjrc, "uint16")
tifffile.imwrite(f"results/tcj{filename}.tif", tcjrc)
tcj = []
for coord in outerTCJ:
tcj.append(findtcjContour(coord, contour[0:-1]))
if "False" in tcj:
tcj.remove("False")
print("removed")
return tcj
def isBoundary(contour):
boundary = False
for con in contour:
if con[0] == 0:
boundary = True
if con[1] == 0:
boundary = True
if con[0] == 511:
boundary = True
if con[1] == 511:
boundary = True
return boundary
def findtcjContour(coord, contour):
close = []
for con in contour:
r = ((con[0] - coord[0]) ** 2 + (con[1] - coord[1]) ** 2) ** 0.5
if r < 1.5:
close.append(con)
if len(close) == 1:
tcj = close[0]
elif len(close) == 0:
tcj = "False"
else:
tcj = np.mean(close, axis=0)
return tcj
def getSecondColour(track, colour):
colours = track[np.all((track - colour) != 0, axis=1)]
colours = colours[np.all((colours - np.array([255, 255, 255])) != 0, axis=1)]
col = []
count = []
while len(colours) > 0:
col.append(colours[0])
count.append(len(colours[np.all((colours - colours[0]) == 0, axis=1)]))
colours = colours[np.all((colours - colours[0]) != 0, axis=1)]
maxm = np.max(count)
colourD = col[count.index(maxm)]
return colourD
def vectorBoundary(pts):
n = len(pts)
test = True
for i in range(n):
if pts[i][0] == 0 or pts[i][0] == 511:
test = False
elif pts[i][1] == 0 or pts[i][1] == 511:
test = False
else:
continue
return test
# -------------------
filenames, fileType = util.getFilesType()
scale = 123.26 / 512
T = 91
# Nucleus velocity relative to tissue
if True:
_df2 = []
_df = []
for filename in filenames:
dist = sm.io.imread(f"dat/{filename}/distance{filename}.tif").astype(int)
df = pd.read_pickle(f"dat/{filename}/shape{filename}.pkl")
Q = np.mean(df["q"])
theta0 = np.arccos(Q[0, 0] / (Q[0, 0] ** 2 + Q[0, 1] ** 2) ** 0.5) / 2
R = util.rotation_matrix(-theta0)
df = pd.read_pickle(f"dat/{filename}/nucleusVelocity{filename}.pkl")
mig = np.zeros(2)
for t in range(90):
dft = df[df["T"] == t]
V = np.mean(dft["Velocity"]) * scale
V = np.matmul(R, V)
_df.append(
{
"Filename": filename,
"T": t,
"V": V,
}
)
dfVelocityMean = pd.DataFrame(_df)
dfVelocityMean.to_pickle(f"databases/dfVelocityMean{fileType}.pkl")
# Cell division relative to tissue and wound
if True:
_df = []
for filename in filenames:
dfDivision = pd.read_pickle(f"dat/{filename}/dfDivision{filename}.pkl")
dfShape = pd.read_pickle(f"dat/{filename}/shape{filename}.pkl")
Q = np.mean(dfShape["q"])
theta0 = 0.5 * np.arctan2(Q[1, 0], Q[0, 0])
t0 = util.findStartTime(filename)
if "Wound" in filename:
dfWound = pd.read_pickle(f"dat/{filename}/woundsite{filename}.pkl")
dist = sm.io.imread(f"dat/{filename}/distance{filename}.tif").astype(int)
for i in range(len(dfDivision)):
t = dfDivision["T"].iloc[i]
(x_w, y_w) = dfWound["Position"].iloc[t]
x = dfDivision["X"].iloc[i]
y = dfDivision["Y"].iloc[i]
ori = (dfDivision["Orientation"].iloc[i] - theta0 * 180 / np.pi) % 180
theta = (np.arctan2(y - y_w, x - x_w) - theta0) * 180 / np.pi
ori_w = (ori - theta) % 180
if ori > 90:
ori = 180 - ori
if ori_w > 90:
ori_w = 180 - ori_w
theta = (np.arctan2(y - y_w, x - x_w) - theta0) * 180 / np.pi
r = dist[t, 512 - y, x]
_df.append(
{
"Filename": filename,
"Label": dfDivision["Label"].iloc[i],
"T": int(t0 + t * 2), # frames are taken every 2 minutes
"X": x * scale,
"Y": y * scale,
"R": r * scale,
"Theta": theta % 360,
"Orientation": ori,
"Orientation Wound": ori_w,
}
)
else:
dfVelocityMean = pd.read_pickle(f"databases/dfVelocityMean{fileType}.pkl")
dfFilename = dfVelocityMean[
dfVelocityMean["Filename"] == filename
].reset_index()
for i in range(len(dfDivision)):
t = dfDivision["T"].iloc[i]
mig = np.sum(
np.stack(np.array(dfFilename.loc[:t, "V"]), axis=0), axis=0
)
xc = 256 * scale + mig[0]
yc = 256 * scale + mig[1]
x = dfDivision["X"].iloc[i] * scale
y = dfDivision["Y"].iloc[i] * scale
r = ((xc - x) ** 2 + (yc - y) ** 2) ** 0.5
ori = (dfDivision["Orientation"].iloc[i] - theta0 * 180 / np.pi) % 180
theta = (np.arctan2(y - yc, x - xc) - theta0) * 180 / np.pi
ori_w = (ori - theta) % 180
if ori_w > 90:
ori_w = 180 - ori_w
if ori > 90:
ori = 180 - ori
_df.append(
{
"Filename": filename,
"Label": dfDivision["Label"].iloc[i],
"T": int(t0 + t * 2), # frames are taken every t2 minutes
"X": x,
"Y": y,
"R": r,
"Theta": theta % 360,
"Orientation": ori,
"Orientation Wound": ori_w,
}
)
dfDivisions = pd.DataFrame(_df)
dfDivisions.to_pickle(f"databases/dfDivisions{fileType}.pkl")
# Cell Shape relative to tissue
if True:
_df2 = []
dfVelocity = pd.read_pickle(f"databases/dfVelocity{fileType}.pkl")
dfVelocityMean = pd.read_pickle(f"databases/dfVelocityMean{fileType}.pkl")
for filename in filenames:
df = pd.read_pickle(f"dat/{filename}/shape{filename}.pkl")
dfFilename = dfVelocityMean[dfVelocityMean["Filename"] == filename]
dist = sm.io.imread(f"dat/{filename}/distance{filename}.tif").astype(int)
mig = np.zeros(2)
Q = np.mean(df["q"])
theta0 = np.arctan2(Q[0, 1], Q[0, 0]) / 2
R = util.rotation_matrix(-theta0)
for t in range(T):
dft = df[df["Time"] == t]
if len(dft) > 0:
Q = np.matmul(R, np.matmul(np.mean(dft["q"]), np.matrix.transpose(R)))
P = np.matmul(R, np.mean(dft["Polar"]))
for i in range(len(dft)):
[x, y] = [
dft["Centroid"].iloc[i][0],
dft["Centroid"].iloc[i][1],
]
r = dist[t, int(512 - y), int(x)]
q = np.matmul(R, np.matmul(dft["q"].iloc[i], np.matrix.transpose(R)))
dq = q - Q
A = dft["Area"].iloc[i] * scale**2
TrQdq = np.trace(np.matmul(Q, dq))
dp = np.matmul(R, dft["Polar"].iloc[i]) - P
[x, y] = np.matmul(R, np.array([x, y]) * scale)
p = np.matmul(R, dft["Polar"].iloc[i])
_df2.append(
{
"Filename": filename,
"T": t,
"X": x - mig[0],
"Y": y - mig[1],
"R": r,
"Centroid": np.array(dft["Centroid"].iloc[i]) * scale,
"dq": dq,
"q": q,
"TrQdq": TrQdq,
"Area": A,
"dp": dp,
"Polar": p,
"Shape Factor": dft["Shape Factor"].iloc[i],
}
)
if t < 90:
mig += np.array(dfFilename["V"][dfFilename["T"] == t])[0]
dfShape = pd.DataFrame(_df2)
dfShape.to_pickle(f"databases/dfShape{fileType}.pkl")
# Cells Divisions and Shape changes
if True:
dfShape = pd.read_pickle(f"databases/dfShape{fileType}.pkl")
_df = []
_dfTrack = []
for filename in filenames:
# print(filename)
df = pd.read_pickle(f"dat/{filename}/dfDivision{filename}.pkl")
tracks = sm.io.imread(f"dat/{filename}/binaryTracks{filename}.tif").astype(int)
binary = sm.io.imread(f"dat/{filename}/binary{filename}.tif").astype(int)
T, X, Y, C = tracks.shape
binary = util.vidrcxy(binary)
tracks = util.vidrcxyRGB(tracks)
tracksDivisions = tracks
for i in range(len(df)):
label = df["Label"].iloc[i]
# if label == 41:
# print(0)
ori = df["Orientation"].iloc[i] % 180
tm = t = int(df["T"].iloc[i])
x = df["X"].iloc[i]
y = df["Y"].iloc[i]
colour = tracks[t, int(x), int(y)]
if np.all((colour - np.array([255, 255, 255])) == 0):
continue
track = tracks[t][np.all((tracks[t] - colour) == 0, axis=2)]
if len(track) > 1500:
continue
finished = False
t_i = t
while finished == False:
t_i += 1
A0 = len(track)
track = tracks[t_i][np.all((tracks[t_i - 1] - colour) == 0, axis=2)]
A1 = len(track[np.all((track - colour) == 0, axis=1)])
if A1 / A0 < 0.65:
finished = True
try:
colourD = getSecondColour(track, colour)
except:
continue
if A1 > 1500:
finished = True
if t_i == T - 1:
finished = True
colourD = getSecondColour(track, colour)
if A1 > 1500:
continue
tc = t_i - 1
if tc - tm < 3:
if tc > 30:
time = np.linspace(tc, tc - 30, 31)
else:
time = np.linspace(tc, 0, tc + 1)
polyList = []
t = int(time[0])
contour = sm.measure.find_contours(
np.all((tracks[t] - colour) == 0, axis=2), level=0
)[0]
poly = sm.measure.approximate_polygon(contour, tolerance=1)
try:
polygon = Polygon(poly)
pts = list(polygon.exterior.coords)
if vectorBoundary(pts):
tcj = findtcj(polygon, binary[t])
else:
tcj = False
theta = cell.orientation_tcj(tcj)
q = cell.qTensor(polygon)
q_tcj = cell.qTensor_tcj(tcj)
_dfTrack.append(
{
"Filename": filename,
"Label": label,
"Type": "parent",
"Colour": colour,
"Time": t,
"Division Time": int(t - tm),
"Polygon": polygon,
"TCJ": tcj,
"Area": polygon.area,
"Shape Factor": cell.shapeFactor(polygon),
"Orientation": (cell.orientation(polygon) * 180 / np.pi)
% 180,
"q": q,
"q0": (q[0, 0] ** 2 + q[0, 1] ** 2) ** 0.5,
"q_tcj": q_tcj,
"q0_tcj": (q_tcj[0, 0] ** 2 + q_tcj[0, 1] ** 2) ** 0.5,
"Shape Factor tcj": cell.shapeFactor_tcj(tcj),
"Orientation tcj": (cell.orientation_tcj(tcj) * 180 / np.pi)
% 180,
}
)
polyList.append(polygon)
polygon0 = polygon
except:
continue
for t in time[1:]:
try:
contour = sm.measure.find_contours(
np.all((tracks[int(t)] - colour) == 0, axis=2), level=0
)[0]
poly = sm.measure.approximate_polygon(contour, tolerance=1)
polygon = Polygon(poly)
if polygon.area == 0:
break
elif dist(polygon, polygon0) > 10:
break
elif polygon0.area / polygon.area < 2 / 3:
break
elif polygon.area / polygon0.area < 2 / 3:
break
elif polygon.area > 1500:
break
else:
pts = list(polygon.exterior.coords)
if vectorBoundary(pts):
tcj = findtcj(polygon, binary[int(t)])
else:
tcj = False
q = cell.qTensor(polygon)
q_tcj = cell.qTensor_tcj(tcj)
_dfTrack.append(
{
"Filename": filename,
"Label": label,
"Type": "parent",
"Colour": colour,
"Time": int(t),
"Division Time": int(t - tm),
"Polygon": polygon,
"TCJ": tcj,
"Area": polygon.area,
"Shape Factor": cell.shapeFactor(polygon),
"Orientation": (
cell.orientation(polygon) * 180 / np.pi
)
% 180,
"q": q,
"q0": (q[0, 0] ** 2 + q[0, 1] ** 2) ** 0.5,
"q_tcj": q_tcj,
"q0_tcj": (q_tcj[0, 0] ** 2 + q_tcj[0, 1] ** 2)
** 0.5,
"Shape Factor tcj": cell.shapeFactor_tcj(tcj),
"Orientation tcj": (
cell.orientation_tcj(tcj) * 180 / np.pi
)
% 180,
}
)
polyList.append(polygon)
polygon0 = polygon
except:
break
if tc < T - 31:
timeD = np.linspace(tc + 1, tc + 31, 31)
else:
timeD = np.linspace(tc + 1, T - 1, T - tc - 1)
polyListD1 = []
t = int(timeD[0])
try:
contour = sm.measure.find_contours(
np.all((tracks[t] - colour) == 0, axis=2), level=0
)[0]
poly = sm.measure.approximate_polygon(contour, tolerance=1)
polygon = Polygon(poly)
if polygon.area == 0:
continue
elif polygon.area > 1500:
continue
else:
pts = list(polygon.exterior.coords)
if vectorBoundary(pts):
tcj = findtcj(polygon, binary[t])
else:
tcj = False
theta = cell.orientation_tcj(tcj)
q = cell.qTensor(polygon)
q_tcj = cell.qTensor_tcj(tcj)
_dfTrack.append(
{
"Filename": filename,
"Label": label,
"Type": "daughter1",
"Colour": colour,
"Time": t,
"Division Time": int(t - tm),
"Polygon": polygon,
"TCJ": tcj,
"Area": polygon.area,
"Shape Factor": cell.shapeFactor(polygon),
"Orientation": (cell.orientation(polygon) * 180 / np.pi)
% 180,
"q": q,
"q0": (q[0, 0] ** 2 + q[0, 1] ** 2) ** 0.5,
"q_tcj": q_tcj,
"q0_tcj": (q_tcj[0, 0] ** 2 + q_tcj[0, 1] ** 2) ** 0.5,
"Shape Factor tcj": cell.shapeFactor_tcj(tcj),
"Orientation tcj": (
cell.orientation_tcj(tcj) * 180 / np.pi
)
% 180,
}
)
polyListD1.append(polygon)
polygon0 = polygon
except:
continue
for t in timeD[1:]:
try:
contour = sm.measure.find_contours(
np.all((tracks[int(t)] - colour) == 0, axis=2), level=0
)[0]
poly = sm.measure.approximate_polygon(contour, tolerance=1)
polygon = Polygon(poly)
if polygon.area == 0:
break
elif dist(polygon, polygon0) > 10:
break
elif polygon0.area / polygon.area < 2 / 3:
break
elif polygon.area / polygon0.area < 2 / 3:
break
elif polygon.area > 1500:
break
else:
pts = list(polygon.exterior.coords)
if vectorBoundary(pts):
tcj = findtcj(polygon, binary[int(t)])
else:
tcj = False
q = cell.qTensor(polygon)
q_tcj = cell.qTensor_tcj(tcj)
_dfTrack.append(
{
"Filename": filename,
"Label": label,
"Type": "daughter1",
"Colour": colour,
"Time": int(t),
"Division Time": int(t - tm),
"Polygon": polygon,
"TCJ": tcj,
"Area": polygon.area,
"Shape Factor": cell.shapeFactor(polygon),
"Orientation": (
cell.orientation(polygon) * 180 / np.pi
)
% 180,
"q": q,
"q0": (q[0, 0] ** 2 + q[0, 1] ** 2) ** 0.5,
"q_tcj": q_tcj,
"q0_tcj": (q_tcj[0, 0] ** 2 + q_tcj[0, 1] ** 2)
** 0.5,
"Shape Factor tcj": cell.shapeFactor_tcj(tcj),
"Orientation tcj": (
cell.orientation_tcj(tcj) * 180 / np.pi
)
% 180,
}
)
polyListD1.append(polygon)
polygon0 = polygon
except:
break
polyListD2 = []
t = int(timeD[0])
try:
contour = sm.measure.find_contours(
np.all((tracks[t] - colourD) == 0, axis=2), level=0
)[0]
poly = sm.measure.approximate_polygon(contour, tolerance=1)
polygon = Polygon(poly)
if polygon.area == 0:
continue
elif polygon.area > 1500:
continue
else:
pts = list(polygon.exterior.coords)
if vectorBoundary(pts):
tcj = findtcj(polygon, binary[t])
else:
tcj = False
theta = cell.orientation_tcj(tcj)
q = cell.qTensor(polygon)
q_tcj = cell.qTensor_tcj(tcj)
_dfTrack.append(
{
"Filename": filename,
"Label": label,
"Type": "daughter2",
"Colour": colourD,
"Time": t,
"Division Time": int(t - tm),
"Polygon": polygon,
"TCJ": tcj,
"Area": polygon.area,
"Shape Factor": cell.shapeFactor(polygon),
"Orientation": (cell.orientation(polygon) * 180 / np.pi)
% 180,
"q": q,
"q0": (q[0, 0] ** 2 + q[0, 1] ** 2) ** 0.5,
"q_tcj": q_tcj,
"q0_tcj": (q_tcj[0, 0] ** 2 + q_tcj[0, 1] ** 2) ** 0.5,
"Shape Factor tcj": cell.shapeFactor_tcj(tcj),
"Orientation tcj": (
cell.orientation_tcj(tcj) * 180 / np.pi
)
% 180,
}
)
polyListD2.append(polygon)
polygon0 = polygon
except:
continue
for t in timeD[1:]:
try:
contour = sm.measure.find_contours(
np.all((tracks[int(t)] - colourD) == 0, axis=2), level=0
)[0]
poly = sm.measure.approximate_polygon(contour, tolerance=1)
polygon = Polygon(poly)
if polygon.area == 0:
break
elif dist(polygon, polygon0) > 10:
break
elif polygon0.area / polygon.area < 2 / 3:
break
elif polygon.area / polygon0.area < 2 / 3:
break
elif polygon.area > 1500:
break
else:
pts = list(polygon.exterior.coords)
if vectorBoundary(pts):
tcj = findtcj(polygon, binary[int(t)])
else:
tcj = False
q = cell.qTensor(polygon)
q_tcj = cell.qTensor_tcj(tcj)
_dfTrack.append(
{
"Filename": filename,
"Label": label,
"Type": "daughter2",
"Colour": colourD,
"Time": int(t),
"Division Time": int(t - tm),
"Polygon": polygon,
"TCJ": tcj,
"Area": polygon.area,
"Shape Factor": cell.shapeFactor(polygon),
"Orientation": (
cell.orientation(polygon) * 180 / np.pi
)
% 180,
"q": q,
"q0": (q[0, 0] ** 2 + q[0, 1] ** 2) ** 0.5,
"q_tcj": q_tcj,
"q0_tcj": (q_tcj[0, 0] ** 2 + q_tcj[0, 1] ** 2)
** 0.5,
"Shape Factor tcj": cell.shapeFactor_tcj(tcj),
"Orientation tcj": (
cell.orientation_tcj(tcj) * 180 / np.pi
)
% 180,
}
)
polyListD2.append(polygon)
polygon0 = polygon
except:
break
_df.append(
{
"Filename": filename,
"Label": label,
"Orientation": ori,
"Times": np.array(time),
"Times daughters": np.array(timeD),
"Cytokineses Time": tc,
"Anaphase Time": tm,
"X": x,
"Y": y,
"Colour": colour,
"Daughter Colour": colourD,
"Polygons": polyList,
"Polygons daughter1": polyListD1,
"Polygons daughter2": polyListD2,
"Track length": len(polyList),
"Daughter length": np.min([len(polyListD1), len(polyListD2)]),
}
)
dfDivisionShape = pd.DataFrame(_df)
dfDivisionShape.to_pickle(f"databases/dfDivisionShape{fileType}.pkl")
dfDivisionTrack = pd.DataFrame(_dfTrack)
dfDivisionTrack.to_pickle(f"databases/dfDivisionTrack{fileType}.pkl")
