JFoz commited on
Commit
016c3b8
1 Parent(s): 305e5eb

Add checks for removed foci

Browse files
path_analysis/analyse.py CHANGED
@@ -36,22 +36,25 @@ def get_paths_from_traces_file(traces_file):
36
  path_lengths.append(float(length))
37
  return all_paths, path_lengths
38
 
39
- def calculate_path_length(point_list, voxel_size=(1,1,1)):
40
- # Simple calculation
41
- l = 0
42
- s = np.array(voxel_size)
43
- for i in range(len(point_list)-1):
44
- l += la.norm(s * (np.array(point_list[i+1]) - np.array(point_list[i])))
45
- return l
46
 
47
 
48
  def calculate_path_length_partials(point_list, voxel_size=(1,1,1)):
 
 
 
 
 
 
 
 
 
 
49
  # Simple calculation
50
- l = [0.0]
51
  s = np.array(voxel_size)
52
  for i in range(len(point_list)-1):
53
- l.append(la.norm(s * (np.array(point_list[i+1]) - np.array(point_list[i]))))
54
- return np.cumsum(l)
55
 
56
 
57
  def visualise_ordering(points_list, dim, wr=5, wc=5):
@@ -260,6 +263,7 @@ def measure_chrom2(path, hei10, config):
260
  measurements = measure_all_with_sphere(path, hei10, op='mean', R=sphere_xy_radius, z_scale_ratio=scale_ratio)
261
  measurements_max = measure_all_with_sphere(path, hei10, op='max', R=sphere_xy_radius, z_scale_ratio=scale_ratio)
262
 
 
263
  return vis, measurements, measurements_max
264
 
265
  def extract_peaks(cell_id, all_paths, path_lengths, measured_traces, config):
@@ -286,7 +290,7 @@ def extract_peaks(cell_id, all_paths, path_lengths, measured_traces, config):
286
  n_paths = len(all_paths)
287
 
288
  data = []
289
- foci_absolute_intensity, foci_position, foci_position_index, trace_median_intensities, trace_thresholds = analyse_traces(all_paths, path_lengths, measured_traces, config)
290
 
291
  foci_intensities = []
292
  for path_foci_abs_int, tmi in zip(foci_absolute_intensity, trace_median_intensities):
@@ -299,6 +303,8 @@ def extract_peaks(cell_id, all_paths, path_lengths, measured_traces, config):
299
 
300
  pl = calculate_path_length_partials(all_paths[i], (config['xy_res'], config['xy_res'], config['z_res']))
301
 
 
 
302
  path_data = { 'Cell_ID':cell_id,
303
  'Trace': i+1,
304
  'SNT_trace_length(um)': path_lengths[i],
@@ -315,7 +321,7 @@ def extract_peaks(cell_id, all_paths, path_lengths, measured_traces, config):
315
  path_data[f'Foci_{j+1}_relative_intensity'] = (v - trace_median_intensities[i])/mean_intensity
316
  data.append(path_data)
317
  trace_positions.append(pl)
318
- return pd.DataFrame(data), foci_absolute_intensity, foci_position_index, trace_thresholds, trace_positions
319
 
320
 
321
  def analyse_paths(cell_id,
@@ -351,25 +357,32 @@ def analyse_paths(cell_id,
351
  vis, m, _ = measure_chrom2(p,foci_stack.transpose(2,1,0), config)
352
  all_trace_vis.append(vis)
353
  all_m.append(m)
 
354
 
355
-
356
- extracted_peaks, foci_absolute_intensity, foci_pos_index, trace_thresholds, trace_positions = extract_peaks(cell_id, all_paths, path_lengths, all_m, config)
357
 
358
 
359
  n_cols = 2
360
  n_rows = (len(all_paths)+n_cols-1)//n_cols
361
- fig, ax = plt.subplots(n_rows,n_cols)
362
  ax = ax.flatten()
363
 
364
  for i, m in enumerate(all_m):
365
  ax[i].set_title(f'Trace {i+1}')
366
  ax[i].plot(trace_positions[i], m)
367
- print(foci_pos_index)
368
  if len(foci_pos_index[i]):
369
  ax[i].plot(trace_positions[i][foci_pos_index[i]], np.array(m)[foci_pos_index[i]], 'rx')
370
- ax[i].set_xlabel('Distance from start (um)')
371
- ax[i].set_ylabel('Intensity')
 
 
 
 
 
 
372
  ax[i].axhline(trace_thresholds[i], c='r', ls=':')
 
 
373
  for i in range(len(all_m), n_cols*n_rows):
374
  ax[i].axis('off')
375
 
 
36
  path_lengths.append(float(length))
37
  return all_paths, path_lengths
38
 
 
 
 
 
 
 
 
39
 
40
 
41
  def calculate_path_length_partials(point_list, voxel_size=(1,1,1)):
42
+ """
43
+ Calculate the partial path length of a series of points.
44
+
45
+ Args:
46
+ point_list (list of tuple): List of points, each represented as a tuple of coordinates (x, y, z).
47
+ voxel_size (tuple, optional): Size of the voxel in each dimension (x, y, z). Defaults to (1, 1, 1).
48
+
49
+ Returns:
50
+ numpy.ndarray: Array of cumulative partial path lengths at each point.
51
+ """
52
  # Simple calculation
53
+ section_lengths = [0.0]
54
  s = np.array(voxel_size)
55
  for i in range(len(point_list)-1):
56
+ section_lengths.append(la.norm(s * (np.array(point_list[i+1]) - np.array(point_list[i]))))
57
+ return np.cumsum(section_lengths)
58
 
59
 
60
  def visualise_ordering(points_list, dim, wr=5, wc=5):
 
263
  measurements = measure_all_with_sphere(path, hei10, op='mean', R=sphere_xy_radius, z_scale_ratio=scale_ratio)
264
  measurements_max = measure_all_with_sphere(path, hei10, op='max', R=sphere_xy_radius, z_scale_ratio=scale_ratio)
265
 
266
+
267
  return vis, measurements, measurements_max
268
 
269
  def extract_peaks(cell_id, all_paths, path_lengths, measured_traces, config):
 
290
  n_paths = len(all_paths)
291
 
292
  data = []
293
+ foci_absolute_intensity, foci_position, foci_position_index, dominated_foci_data, trace_median_intensities, trace_thresholds = analyse_traces(all_paths, path_lengths, measured_traces, config)
294
 
295
  foci_intensities = []
296
  for path_foci_abs_int, tmi in zip(foci_absolute_intensity, trace_median_intensities):
 
303
 
304
  pl = calculate_path_length_partials(all_paths[i], (config['xy_res'], config['xy_res'], config['z_res']))
305
 
306
+ print(i, len(all_paths[i]), len(pl))
307
+
308
  path_data = { 'Cell_ID':cell_id,
309
  'Trace': i+1,
310
  'SNT_trace_length(um)': path_lengths[i],
 
321
  path_data[f'Foci_{j+1}_relative_intensity'] = (v - trace_median_intensities[i])/mean_intensity
322
  data.append(path_data)
323
  trace_positions.append(pl)
324
+ return pd.DataFrame(data), foci_absolute_intensity, foci_position_index, dominated_foci_data, trace_thresholds, trace_positions
325
 
326
 
327
  def analyse_paths(cell_id,
 
357
  vis, m, _ = measure_chrom2(p,foci_stack.transpose(2,1,0), config)
358
  all_trace_vis.append(vis)
359
  all_m.append(m)
360
+
361
 
362
+ extracted_peaks, foci_absolute_intensity, foci_pos_index, dominated_foci_data, trace_thresholds, trace_positions = extract_peaks(cell_id, all_paths, path_lengths, all_m, config)
 
363
 
364
 
365
  n_cols = 2
366
  n_rows = (len(all_paths)+n_cols-1)//n_cols
367
+ fig, ax = plt.subplots(n_rows,n_cols, figsize=(5*n_cols, 3*n_rows))
368
  ax = ax.flatten()
369
 
370
  for i, m in enumerate(all_m):
371
  ax[i].set_title(f'Trace {i+1}')
372
  ax[i].plot(trace_positions[i], m)
 
373
  if len(foci_pos_index[i]):
374
  ax[i].plot(trace_positions[i][foci_pos_index[i]], np.array(m)[foci_pos_index[i]], 'rx')
375
+
376
+ if len(dominated_foci_data[i]):
377
+
378
+ dominated_foci_pos_index = [u.idx for u in dominated_foci_data[i]]
379
+ ax[i].plot(trace_positions[i][dominated_foci_pos_index], np.array(m)[dominated_foci_pos_index], color=(0.5,0.5,0.5), marker='o', linestyle='None')
380
+
381
+
382
+ if trace_thresholds[i] is not None:
383
  ax[i].axhline(trace_thresholds[i], c='r', ls=':')
384
+ ax[i].set_xlabel('Distance from start (um)')
385
+ ax[i].set_ylabel('Intensity')
386
  for i in range(len(all_m), n_cols*n_rows):
387
  ax[i].axis('off')
388
 
path_analysis/data_preprocess.py CHANGED
@@ -9,45 +9,54 @@ from math import ceil
9
 
10
 
11
 
12
- def thin_points(point_list, dmin=10, voxel_size=(1,1,1)):
13
  """
14
- Remove points within a specified distance of each other, retaining the point with the highest intensity.
15
 
16
  Args:
17
- - point_list (list of tuples): Each tuple contains:
18
- - x (list of float): 3D coordinates of the point.
19
- - intensity (float): The intensity value of the point.
20
- - idx (int): A unique identifier or index for the point.
21
- - dmin (float, optional): Minimum distance between points. Points closer than this threshold will be thinned. Defaults to 10.
 
22
 
23
  Returns:
24
- - list of int: A list containing indices of the removed points.
 
 
25
 
26
  Notes:
27
- - The function uses the L2 norm (Euclidean distance) to compute the distance between points.
28
- - When two points are within `dmin` distance, the point with the lower intensity is removed.
29
  """
30
- removed_points = []
31
- for i in range(len(point_list)):
32
- if point_list[i][2] in removed_points:
 
33
  continue
34
- for j in range(len(point_list)):
35
  if i==j:
36
  continue
37
- if point_list[j][2] in removed_points:
38
  continue
39
- d = (np.array(point_list[i][0]) - np.array(point_list[j][0]))*np.array(voxel_size)
40
  d = la.norm(d)
41
  if d<dmin:
42
- hi = point_list[i][1]
43
- hj = point_list[j][1]
44
  if hi<hj:
45
- removed_points.append(point_list[i][2])
 
46
  break
47
  else:
48
- removed_points.append(point_list[j][2])
49
-
50
- return removed_points
 
 
 
 
51
 
52
 
53
  @dataclass
@@ -59,6 +68,17 @@ class CellData(object):
59
  """
60
  pathdata_list: list
61
 
 
 
 
 
 
 
 
 
 
 
 
62
  @dataclass
63
  class PathData(object):
64
  """Represents data related to a specific path in the cell.
@@ -67,16 +87,23 @@ class PathData(object):
67
  the defining points, the fluorescence values, and the path length of a specific path.
68
 
69
  Attributes: peaks (list): List of peaks in the path (indicies of positions in points, o_hei10).
 
70
  points (list): List of points defining the path.
71
  o_hei10 (list): List of (unnormalized) fluorescence intensity values along the path
72
  SC_length (float): Length of the path.
73
 
74
  """
75
  peaks: list
 
76
  points: list
77
  o_hei10: list
78
  SC_length: float
79
 
 
 
 
 
 
80
 
81
 
82
  def find_peaks2(v, distance=5, prominence=0.5):
@@ -136,7 +163,7 @@ def process_cell_traces(all_paths, path_lengths, measured_trace_fluorescence):
136
 
137
  cell_peaks = []
138
 
139
- for points, path_length, o_hei10 in zip(all_paths, path_lengths, measured_trace_fluorescence):
140
 
141
  # For peak determination normalize each trace to have mean zero and s.d. 1
142
  hei10_normalized = (o_hei10 - np.mean(o_hei10))/np.std(o_hei10)
@@ -158,24 +185,34 @@ def process_cell_traces(all_paths, path_lengths, measured_trace_fluorescence):
158
  to_thin = []
159
  for k in range(len(cell_peaks)):
160
  for u in range(len(cell_peaks[k][0])):
161
- to_thin.append((cell_peaks[k][1][u], cell_peaks[k][2][u], (k, u)))
162
 
163
  # Exclude any peak with a nearby brighter peak (on any SC)
164
- removed_points = thin_points(to_thin)
165
 
166
-
167
  # Clean up and remove these peaks
168
  new_cell_peaks = []
169
- for k in range(len(cell_peaks)):
170
- cc = []
171
- pp = cell_peaks[k][0]
172
- for u in range(len(pp)):
173
- if (k,u) not in removed_points:
174
- cc.append(pp[u])
175
- new_cell_peaks.append(cc)
 
 
 
 
 
 
 
 
 
 
 
 
176
 
177
  cell_peaks = new_cell_peaks
178
-
179
  pd_list = []
180
 
181
  # Save peak positions, absolute HEI10 intensities, and length for each SC
@@ -184,8 +221,9 @@ def process_cell_traces(all_paths, path_lengths, measured_trace_fluorescence):
184
  points, o_hei10 = all_paths[k], measured_trace_fluorescence[k]
185
 
186
  peaks = cell_peaks[k]
 
187
 
188
- pd = PathData(peaks=peaks, points=points, o_hei10=o_hei10, SC_length=path_lengths[k])
189
  pd_list.append(pd)
190
 
191
  cd = CellData(pathdata_list=pd_list)
@@ -196,9 +234,9 @@ def process_cell_traces(all_paths, path_lengths, measured_trace_fluorescence):
196
  alpha_max = 0.4
197
 
198
 
199
- # Criterion used for identifying peak as a CO - normalized (with mean and s.d.)
200
  # hei10 levels being above 0.4 time maximum peak level
201
- def pc(pos, v, alpha=alpha_max):
202
  """
203
  Identify and return positions where values in the array `v` exceed a certain threshold.
204
 
@@ -212,8 +250,11 @@ def pc(pos, v, alpha=alpha_max):
212
  Returns:
213
  - numpy.ndarray: Array of positions where corresponding values in `v` exceed the threshold.
214
  """
215
- idx = (v>=alpha*np.max(v))
216
- return np.array(pos[idx])
 
 
 
217
 
218
  def analyse_celldata(cell_data, config):
219
  """
@@ -226,14 +267,18 @@ def analyse_celldata(cell_data, config):
226
  'threshold_type' (str) = 'per-trace', 'per-foci'
227
 
228
  Returns:
229
- tuple: A tuple containing three lists:
230
  - foci_rel_intensity (list): List of relative intensities for the detected foci.
231
  - foci_pos (list): List of absolute positions of the detected foci.
232
  - foci_pos_index (list): List of indices of the detected foci.
 
 
 
233
  """
234
  foci_abs_intensity = []
235
  foci_pos = []
236
  foci_pos_index = []
 
237
  trace_median_intensities = []
238
  trace_thresholds = []
239
 
@@ -254,15 +299,39 @@ def analyse_celldata(cell_data, config):
254
  h = np.array(path_data.o_hei10)
255
  h = h - np.mean(h)
256
  h = h/np.std(h)
257
- # Extract peaks according to criterion
258
- sig_peak_idx = pc(peaks, h[peaks], peak_threshold)
259
- trace_thresholds.append((1-peak_threshold)*np.mean(path_data.o_hei10) + peak_threshold*np.max(np.array(path_data.o_hei10)[peaks]))
 
 
 
 
 
260
 
261
- pos_abs = (sig_peak_idx/len(path_data.points))*path_data.SC_length
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
262
  foci_pos.append(pos_abs)
263
- foci_abs_intensity.append(np.array(path_data.o_hei10)[sig_peak_idx])
264
 
265
- foci_pos_index.append(sig_peak_idx)
266
  trace_median_intensities.append(np.median(path_data.o_hei10))
267
 
268
  elif threshold_type == 'per-cell':
@@ -286,22 +355,34 @@ def analyse_celldata(cell_data, config):
286
  h = np.array(path_data.o_hei10)
287
  h = h - np.mean(h)
288
 
289
- sig_peak_idx = peaks[h[peaks]>peak_threshold*max_cell_intensity]
 
 
 
290
 
291
  trace_thresholds.append(np.mean(path_data.o_hei10) + peak_threshold*max_cell_intensity)
292
 
 
 
 
 
 
 
 
 
 
293
 
294
- pos_abs = (sig_peak_idx/len(path_data.points))*path_data.SC_length
295
  foci_pos.append(pos_abs)
296
- foci_abs_intensity.append(np.array(path_data.o_hei10)[sig_peak_idx])
297
 
298
- foci_pos_index.append(sig_peak_idx)
299
  trace_median_intensities.append(np.median(path_data.o_hei10))
300
 
301
  else:
302
  raise NotImplementedError
303
 
304
- return foci_abs_intensity, foci_pos, foci_pos_index, trace_median_intensities, trace_thresholds
305
 
306
  def analyse_traces(all_paths, path_lengths, measured_trace_fluorescence, config):
307
 
 
9
 
10
 
11
 
12
+ def thin_peaks(peak_list, dmin=10, voxel_size=(1,1,1), return_larger_peaks=False):
13
  """
14
+ Remove peaks within a specified distance of each other, retaining the peak with the highest intensity.
15
 
16
  Args:
17
+ - peak_list (list of PeakData): Each element contains:
18
+ - pos (list of float): 3D coordinates of the peak.
19
+ - intensity (float): The intensity value of the peak.
20
+ - key (tuple): A unique identifier or index for the peak (#trace, #peak)
21
+ - dmin (float, optional): Minimum distance between peaks. peaks closer than this threshold will be thinned. Defaults to 10.
22
+ - return_larger_peaks (bool, optional): Indicate larger peak for each thinned peak
23
 
24
  Returns:
25
+ - list of tuples: A list containing keys of the removed peaks.
26
+ if return_larger_peaks
27
+ - list of tuples: A list containing the keys of the larger peak causing the peak to be removed
28
 
29
  Notes:
30
+ - The function uses the L2 norm (Euclidean distance) to compute the distance between peaks.
31
+ - When two peaks are within `dmin` distance, the peak with the lower intensity is removed.
32
  """
33
+ removed_peaks = []
34
+ removed_larger_peaks = []
35
+ for i in range(len(peak_list)):
36
+ if peak_list[i].key in removed_peaks:
37
  continue
38
+ for j in range(len(peak_list)):
39
  if i==j:
40
  continue
41
+ if peak_list[j].key in removed_peaks:
42
  continue
43
+ d = (np.array(peak_list[i].pos) - np.array(peak_list[j].pos))*np.array(voxel_size)
44
  d = la.norm(d)
45
  if d<dmin:
46
+ hi = peak_list[i].intensity
47
+ hj = peak_list[j].intensity
48
  if hi<hj:
49
+ removed_peaks.append(peak_list[i].key)
50
+ removed_larger_peaks.append(peak_list[j].key)
51
  break
52
  else:
53
+ removed_peaks.append(peak_list[j].key)
54
+ removed_larger_peaks.append(peak_list[i].key)
55
+
56
+ if return_larger_peaks:
57
+ return removed_peaks, removed_larger_peaks
58
+ else:
59
+ return removed_peaks
60
 
61
 
62
  @dataclass
 
68
  """
69
  pathdata_list: list
70
 
71
+ @dataclass
72
+ class RemovedPeakData(object):
73
+ """Represents data related to a removed peak
74
+
75
+ Attributes:
76
+ idx (int): Index of peak along path
77
+ dominating_peak (tuple): (path_idx, position along path) for dominating peak
78
+ """
79
+ idx: int
80
+ dominating_peak: tuple
81
+
82
  @dataclass
83
  class PathData(object):
84
  """Represents data related to a specific path in the cell.
 
87
  the defining points, the fluorescence values, and the path length of a specific path.
88
 
89
  Attributes: peaks (list): List of peaks in the path (indicies of positions in points, o_hei10).
90
+ removed_peaks (list): List of peaks in the path which have been removed because of a nearby larger peak
91
  points (list): List of points defining the path.
92
  o_hei10 (list): List of (unnormalized) fluorescence intensity values along the path
93
  SC_length (float): Length of the path.
94
 
95
  """
96
  peaks: list
97
+ removed_peaks: list
98
  points: list
99
  o_hei10: list
100
  SC_length: float
101
 
102
+ @dataclass
103
+ class PeakData(object):
104
+ pos: tuple
105
+ intensity: float
106
+ key: tuple
107
 
108
 
109
  def find_peaks2(v, distance=5, prominence=0.5):
 
163
 
164
  cell_peaks = []
165
 
166
+ for points, o_hei10 in zip(all_paths, measured_trace_fluorescence):
167
 
168
  # For peak determination normalize each trace to have mean zero and s.d. 1
169
  hei10_normalized = (o_hei10 - np.mean(o_hei10))/np.std(o_hei10)
 
185
  to_thin = []
186
  for k in range(len(cell_peaks)):
187
  for u in range(len(cell_peaks[k][0])):
188
+ to_thin.append(PeakData(pos=cell_peaks[k][1][u], intensity=cell_peaks[k][2][u], key=(k, u)))
189
 
190
  # Exclude any peak with a nearby brighter peak (on any SC)
191
+ removed_peaks, removed_larger_peaks = thin_peaks(to_thin, return_larger_peaks=True)
192
 
 
193
  # Clean up and remove these peaks
194
  new_cell_peaks = []
195
+ removed_cell_peaks = []
196
+ removed_cell_peaks_larger = []
197
+ for path_idx in range(len(cell_peaks)):
198
+ path_retained_peaks = []
199
+ path_removed_peaks = []
200
+ path_peaks = cell_peaks[path_idx][0]
201
+
202
+ for peak_idx in range(len(path_peaks)):
203
+ if (path_idx, peak_idx) not in removed_peaks:
204
+ path_retained_peaks.append(path_peaks[peak_idx])
205
+ else:
206
+ # What's the larger point?
207
+ idx = removed_peaks.index((path_idx, peak_idx))
208
+ larger_path, larger_idx = removed_larger_peaks[idx]
209
+ path_removed_peaks.append(RemovedPeakData(idx=path_peaks[peak_idx], dominating_peak=(larger_path, cell_peaks[larger_path][0][larger_idx])))
210
+ ###
211
+
212
+ new_cell_peaks.append(path_retained_peaks)
213
+ removed_cell_peaks.append(path_removed_peaks)
214
 
215
  cell_peaks = new_cell_peaks
 
216
  pd_list = []
217
 
218
  # Save peak positions, absolute HEI10 intensities, and length for each SC
 
221
  points, o_hei10 = all_paths[k], measured_trace_fluorescence[k]
222
 
223
  peaks = cell_peaks[k]
224
+ removed_peaks = removed_cell_peaks[k]
225
 
226
+ pd = PathData(peaks=peaks, removed_peaks=removed_peaks, points=points, o_hei10=o_hei10, SC_length=path_lengths[k])
227
  pd_list.append(pd)
228
 
229
  cd = CellData(pathdata_list=pd_list)
 
234
  alpha_max = 0.4
235
 
236
 
237
+ # Criterion used for identifying peak as a focus - normalized (with mean and s.d.)
238
  # hei10 levels being above 0.4 time maximum peak level
239
+ def focus_criterion(pos, v, alpha=alpha_max):
240
  """
241
  Identify and return positions where values in the array `v` exceed a certain threshold.
242
 
 
250
  Returns:
251
  - numpy.ndarray: Array of positions where corresponding values in `v` exceed the threshold.
252
  """
253
+ if len(v):
254
+ idx = (v>=alpha*np.max(v))
255
+ return np.array(pos[idx])
256
+ else:
257
+ return np.array([], dtype=np.int32)
258
 
259
  def analyse_celldata(cell_data, config):
260
  """
 
267
  'threshold_type' (str) = 'per-trace', 'per-foci'
268
 
269
  Returns:
270
+ tuple: A tuple containing:
271
  - foci_rel_intensity (list): List of relative intensities for the detected foci.
272
  - foci_pos (list): List of absolute positions of the detected foci.
273
  - foci_pos_index (list): List of indices of the detected foci.
274
+ - dominated_foci_data (list): List of RemovedPeakData indicating positions of removed peaks and the index of the larger peak
275
+ - trace_median_intensities (list): Per-trace median intensity
276
+ - trace_thresholds (list): Per-trace absolute threshold for calling peaks as foci
277
  """
278
  foci_abs_intensity = []
279
  foci_pos = []
280
  foci_pos_index = []
281
+ dominated_foci_data = []
282
  trace_median_intensities = []
283
  trace_thresholds = []
284
 
 
299
  h = np.array(path_data.o_hei10)
300
  h = h - np.mean(h)
301
  h = h/np.std(h)
302
+ # Extract foci according to criterion
303
+ foci_idx = focus_criterion(peaks, h[peaks], peak_threshold)
304
+ print('peaks', peaks, h[peaks], foci_idx, np.mean(path_data.o_hei10))
305
+
306
+ #
307
+ removed_peaks = path_data.removed_peaks
308
+ removed_peaks_idx = np.array([u.idx for u in removed_peaks], dtype=np.int32)
309
+
310
 
311
+ if len(peaks):
312
+ trace_thresholds.append((1-peak_threshold)*np.mean(path_data.o_hei10) + peak_threshold*np.max(np.array(path_data.o_hei10)[peaks]))
313
+ else:
314
+ trace_thresholds.append(None)
315
+
316
+ if len(removed_peaks):
317
+ if len(peaks):
318
+ threshold = (1-peak_threshold)*np.mean(path_data.o_hei10) + peak_threshold*np.max(np.array(path_data.o_hei10)[peaks])
319
+ else:
320
+ threshold = float('-inf')
321
+
322
+
323
+ removed_peak_heights = np.array(path_data.o_hei10)[removed_peaks_idx]
324
+ dominated_foci_idx = np.where(removed_peak_heights>threshold)[0]
325
+
326
+ dominated_foci_data.append([removed_peaks[i] for i in dominated_foci_idx])
327
+ else:
328
+ dominated_foci_data.append([])
329
+
330
+ pos_abs = (foci_idx/len(path_data.points))*path_data.SC_length
331
  foci_pos.append(pos_abs)
332
+ foci_abs_intensity.append(np.array(path_data.o_hei10)[foci_idx])
333
 
334
+ foci_pos_index.append(foci_idx)
335
  trace_median_intensities.append(np.median(path_data.o_hei10))
336
 
337
  elif threshold_type == 'per-cell':
 
355
  h = np.array(path_data.o_hei10)
356
  h = h - np.mean(h)
357
 
358
+ foci_idx = peaks[h[peaks]>peak_threshold*max_cell_intensity]
359
+
360
+ removed_peaks = path_data.removed_peaks
361
+ removed_peaks_idx = np.array([u.idx for u in removed_peaks], dtype=np.int32)
362
 
363
  trace_thresholds.append(np.mean(path_data.o_hei10) + peak_threshold*max_cell_intensity)
364
 
365
+ if len(removed_peaks):
366
+ threshold = np.mean(path_data.o_hei10) + peak_threshold*max_cell_intensity
367
+
368
+ removed_peak_heights = np.array(path_data.o_hei10)[removed_peaks_idx]
369
+ dominated_foci_idx = np.where(removed_peak_heights>threshold)[0]
370
+
371
+ dominated_foci_data.append([removed_peaks[i] for i in dominated_foci_idx])
372
+ else:
373
+ dominated_foci_data.append([])
374
 
375
+ pos_abs = (foci_idx/len(path_data.points))*path_data.SC_length
376
  foci_pos.append(pos_abs)
377
+ foci_abs_intensity.append(np.array(path_data.o_hei10)[foci_idx])
378
 
379
+ foci_pos_index.append(foci_idx)
380
  trace_median_intensities.append(np.median(path_data.o_hei10))
381
 
382
  else:
383
  raise NotImplementedError
384
 
385
+ return foci_abs_intensity, foci_pos, foci_pos_index, dominated_foci_data, trace_median_intensities, trace_thresholds
386
 
387
  def analyse_traces(all_paths, path_lengths, measured_trace_fluorescence, config):
388
 
tests/test_analyse.py CHANGED
@@ -1,9 +1,68 @@
1
 
 
2
  from path_analysis.analyse import *
3
  import numpy as np
4
  from math import pi
5
  import xml.etree.ElementTree as ET
 
6
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
7
 
8
  def test_get_paths_from_traces_file():
9
  # Mock the XML traces file content
@@ -221,3 +280,34 @@ def test_make_sphere_equal():
221
  assert abs(np.sum(sphere)-4/3*pi*R**3)<10, f"Expected approximate volume to be correct"
222
  assert (sphere[R,R,0] == 1), f"Expected centre point on top plane to be within sphere"
223
  assert (sphere[R+1,R,0] == 0), f"Expected point next to centre on top plane to be outside sphere"
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
 
2
+ import pytest
3
  from path_analysis.analyse import *
4
  import numpy as np
5
  from math import pi
6
  import xml.etree.ElementTree as ET
7
+ from PIL import ImageChops
8
 
9
+ def test_draw_paths_no_error():
10
+ all_paths = [[[0, 0], [1, 1]], [[2, 2], [3, 3]]]
11
+ foci_stack = np.zeros((5, 5, 5))
12
+ foci_stack[0,0,0] = 1.0
13
+ foci_index = [[0], [1]]
14
+ r = 3
15
+
16
+ try:
17
+ im = draw_paths(all_paths, foci_stack, foci_index, r)
18
+ except Exception as e:
19
+ pytest.fail(f"draw_paths raised an exception: {e}")
20
+
21
+ def test_draw_paths_image_size():
22
+ all_paths = [[[0, 0], [1, 1]], [[2, 2], [3, 3]]]
23
+ foci_stack = np.zeros((5, 5, 5))
24
+ foci_stack[0,0,0] = 1.0
25
+
26
+ foci_index = [[0], [1]]
27
+ r = 3
28
+
29
+ im = draw_paths(all_paths, foci_stack, foci_index, r)
30
+ assert im.size == (5, 5), f"Expected image size (5, 5), got {im.size}"
31
+
32
+ def test_draw_paths_image_modified():
33
+ all_paths = [[[0, 0], [1, 1]], [[2, 2], [3, 3]]]
34
+ foci_stack = np.zeros((5, 5, 5))
35
+ foci_stack[0,0,0] = 1.0
36
+ foci_index = [[0], [1]]
37
+ r = 3
38
+
39
+ im = draw_paths(all_paths, foci_stack, foci_index, r)
40
+ blank_image = Image.new("RGB", (5, 5), "black")
41
+
42
+ # Check if the image is not entirely black (i.e., has been modified)
43
+ diff = ImageChops.difference(im, blank_image)
44
+ assert diff.getbbox() is not None, "The image has not been modified"
45
+
46
+
47
+
48
+ def test_calculate_path_length_partials_default_voxel():
49
+ point_list = [(0, 0, 0), (1, 0, 0), (1, 1, 1)]
50
+ expected_result = np.array([0.0, 1.0, 1.0+np.sqrt(2)])
51
+ result = calculate_path_length_partials(point_list)
52
+ np.testing.assert_allclose(result, expected_result, atol=1e-5)
53
+
54
+ def test_calculate_path_length_partials_custom_voxel():
55
+ point_list = [(0, 0, 0), (1, 0, 0), (1, 1, 0)]
56
+ voxel_size = (1, 2, 1)
57
+ expected_result = np.array([0.0, 1.0, 3.0])
58
+ result = calculate_path_length_partials(point_list, voxel_size=voxel_size)
59
+ np.testing.assert_allclose(result, expected_result, atol=1e-5)
60
+
61
+ def test_calculate_path_length_partials_single_point():
62
+ point_list = [(0, 0, 0)]
63
+ expected_result = np.array([0.0])
64
+ result = calculate_path_length_partials(point_list)
65
+ np.testing.assert_allclose(result, expected_result, atol=1e-5)
66
 
67
  def test_get_paths_from_traces_file():
68
  # Mock the XML traces file content
 
280
  assert abs(np.sum(sphere)-4/3*pi*R**3)<10, f"Expected approximate volume to be correct"
281
  assert (sphere[R,R,0] == 1), f"Expected centre point on top plane to be within sphere"
282
  assert (sphere[R+1,R,0] == 0), f"Expected point next to centre on top plane to be outside sphere"
283
+
284
+ import pandas as pd
285
+
286
+
287
+ # 1. Test basic functionality
288
+ def test_extract_peaks_basic():
289
+ cell_id = 1
290
+ all_paths = [[[0, 0], [1, 1]]]
291
+ path_lengths = [1.41] # length of the above path
292
+ measured_traces = [[100, 200]] # fluorescence along the path
293
+ config = {'peak_threshold': 0.4, 'sphere_radius': 2, 'xy_res': 1, 'z_res': 1, 'use_corrected_positions': True}
294
+
295
+ df, foci_abs_int, foci_pos_idx, _, _, _ = extract_peaks(cell_id, all_paths, path_lengths, measured_traces, config)
296
+
297
+ # Now add your assertions to validate the result
298
+ assert len(df) == 1, "Expected one row in DataFrame"
299
+ assert df['Cell_ID'].iloc[0] == cell_id, "Unexpected cell_id"
300
+ # Add more assertions here based on expected values
301
+
302
+ # 2. Test multiple paths
303
+ def test_extract_peaks_multiple_paths():
304
+ cell_id = 1
305
+ all_paths = [[[0, 0], [1, 1]], [[1, 1], [2, 2]]]
306
+ path_lengths = [1.41, 1.41]
307
+ measured_traces = [[100, 200], [100, 150]]
308
+ config = {'peak_threshold': 0.4, 'sphere_radius': 2, 'xy_res': 1, 'z_res': 1, 'use_corrected_positions': True}
309
+
310
+ df, _, _, _, _, _ = extract_peaks(cell_id, all_paths, path_lengths, measured_traces, config)
311
+
312
+ assert len(df) == 2, "Expected two rows in DataFrame"
313
+ # Add more assertions here
tests/test_preprocess.py CHANGED
@@ -2,18 +2,19 @@ from path_analysis.data_preprocess import *
2
  import numpy as np
3
  import pytest
4
 
 
5
  def test_thin_points():
6
  # Define a sample point list
7
  points = [
8
- ([0, 0, 0], 10, 0),
9
- ([1, 1, 1], 8, 1),
10
- ([10, 10, 10], 12, 2),
11
- ([10.5, 10.5, 10.5], 5, 3),
12
- ([20, 20, 20], 15, 4)
13
  ]
14
 
15
  # Call the thin_points function with dmin=5 (for example)
16
- removed_indices = thin_points(points, dmin=5)
17
 
18
  # Check results
19
  # Point at index 1 ([1, 1, 1]) should be removed since it's within 5 units distance of point at index 0 and has lower intensity.
@@ -22,12 +23,12 @@ def test_thin_points():
22
 
23
  # Another simple test to check if function does nothing when points are far apart
24
  far_points = [
25
- ([0, 0, 0], 10, 0),
26
- ([100, 100, 100], 12, 1),
27
- ([200, 200, 200], 15, 2)
28
  ]
29
 
30
- removed_indices_far = thin_points(far_points, dmin=5)
31
  assert len(removed_indices_far) == 0 # Expect no points to be removed
32
 
33
 
@@ -72,29 +73,32 @@ def test_find_peaks2():
72
  assert peaks == [1] # Only the peak at position 1 meets the prominence threshold
73
 
74
 
75
- def test_pc():
76
  pos = np.array([0, 1, 2, 3, 4, 6])
77
  values = np.array([0.1, 0.5, 0.2, 0.8, 0.3, 0.9])
78
 
79
  # Basic test
80
- assert np.array_equal(pc(pos, values), np.array([1, 3, 6])) # only values 0.8 and 0.9 exceed 0.4 times the max (which is 0.9)
 
 
 
81
 
82
  # Test with custom alpha
83
- assert np.array_equal(pc(pos, values, alpha=0.5), np.array([1, 3, 6]))
84
 
85
  # Test with a larger alpha
86
- assert np.array_equal(pc(pos, values, alpha=1.0), [6]) # No values exceed the maximum value itself
87
 
88
  # Test with all values below threshold
89
  values = np.array([0.1, 0.2, 0.3, 0.4])
90
 
91
- assert np.array_equal(pc(pos[:4], values), [1,2,3]) # All values are below 0.4 times the max (which is 0.4)
92
 
93
  @pytest.fixture
94
  def mock_data():
95
  all_paths = [ [ (0,0,0), (0,2,0), (0,5,0), (0,10,0), (0,15,0), (0,20,0)], [ (1,20,0), (1,20,10), (1,20,20) ] ] # Mock paths
96
  path_lengths = [ 2.2, 2.3 ] # Mock path lengths
97
- measured_trace_fluorescence = [ [100, 8, 3, 2, 3, 39], [38, 2, 20] ] # Mock fluorescence data
98
  return all_paths, path_lengths, measured_trace_fluorescence
99
 
100
  def test_process_cell_traces_return_type(mock_data):
@@ -117,23 +121,30 @@ def test_process_cell_traces_pathdata_path_lengths(mock_data):
117
  def test_process_cell_traces_peaks(mock_data):
118
  all_paths, path_lengths, measured_trace_fluorescence = mock_data
119
  result = process_cell_traces(all_paths, path_lengths, measured_trace_fluorescence)
 
120
  peaks = [p.peaks for p in result.pathdata_list]
121
  assert peaks == [[0,5],[]]
122
 
123
  # Mock data
124
  @pytest.fixture
125
  def mock_celldata():
126
- pathdata1 = PathData(peaks=[0, 5], points=[(0,0,0), (0,2,0), (0,5,0), (0,10,0), (0,15,0), (0,20,0)], o_hei10=[100, 8, 3, 2, 3, 39], SC_length=2.2)
127
- pathdata2 = PathData(peaks=[0], points=[(1,20,0), (1,20,10), (1,20,20) ], o_hei10=[38, 2, 20], SC_length=2.3)
128
  return CellData(pathdata_list=[pathdata1, pathdata2])
129
 
130
- def test_analyse_celldata_output_length(mock_celldata):
131
- rel_intensity, pos, pos_index, trace_median_intensity, trace_thresholds = analyse_celldata(mock_celldata, {'peak_threshold': 0.4, 'threshold_type':'per-trace'})
132
- assert len(rel_intensity) == len(mock_celldata.pathdata_list), "Mismatch in relative intensities length"
133
- assert len(pos) == len(mock_celldata.pathdata_list), "Mismatch in positions length"
134
- assert len(pos_index) == len(mock_celldata.pathdata_list), "Mismatch in position indices length"
135
-
136
 
 
137
 
138
 
 
 
 
 
 
 
139
 
 
2
  import numpy as np
3
  import pytest
4
 
5
+
6
  def test_thin_points():
7
  # Define a sample point list
8
  points = [
9
+ PeakData([0, 0, 0], 10, 0),
10
+ PeakData([1, 1, 1], 8, 1),
11
+ PeakData([10, 10, 10], 12, 2),
12
+ PeakData([10.5, 10.5, 10.5], 5, 3),
13
+ PeakData([20, 20, 20], 15, 4)
14
  ]
15
 
16
  # Call the thin_points function with dmin=5 (for example)
17
+ removed_indices = thin_peaks(points, dmin=5)
18
 
19
  # Check results
20
  # Point at index 1 ([1, 1, 1]) should be removed since it's within 5 units distance of point at index 0 and has lower intensity.
 
23
 
24
  # Another simple test to check if function does nothing when points are far apart
25
  far_points = [
26
+ PeakData([0, 0, 0], 10, 0),
27
+ PeakData([100, 100, 100], 12, 1),
28
+ PeakData([200, 200, 200], 15, 2)
29
  ]
30
 
31
+ removed_indices_far = thin_peaks(far_points, dmin=5)
32
  assert len(removed_indices_far) == 0 # Expect no points to be removed
33
 
34
 
 
73
  assert peaks == [1] # Only the peak at position 1 meets the prominence threshold
74
 
75
 
76
+ def test_focus_criterion():
77
  pos = np.array([0, 1, 2, 3, 4, 6])
78
  values = np.array([0.1, 0.5, 0.2, 0.8, 0.3, 0.9])
79
 
80
  # Basic test
81
+ assert np.array_equal(focus_criterion(pos, values), np.array([1, 3, 6])) # only values 0.8 and 0.9 exceed 0.4 times the max (which is 0.9)
82
+
83
+ # Empty test
84
+ assert np.array_equal(focus_criterion(np.array([]), np.array([])), np.array([]))
85
 
86
  # Test with custom alpha
87
+ assert np.array_equal(focus_criterion(pos, values, alpha=0.5), np.array([1, 3, 6]))
88
 
89
  # Test with a larger alpha
90
+ assert np.array_equal(focus_criterion(pos, values, alpha=1.0), [6]) # No values exceed the maximum value itself
91
 
92
  # Test with all values below threshold
93
  values = np.array([0.1, 0.2, 0.3, 0.4])
94
 
95
+ assert np.array_equal(focus_criterion(pos[:4], values), [1,2,3]) # All values are below 0.4 times the max (which is 0.4)
96
 
97
  @pytest.fixture
98
  def mock_data():
99
  all_paths = [ [ (0,0,0), (0,2,0), (0,5,0), (0,10,0), (0,15,0), (0,20,0)], [ (1,20,0), (1,20,10), (1,20,20) ] ] # Mock paths
100
  path_lengths = [ 2.2, 2.3 ] # Mock path lengths
101
+ measured_trace_fluorescence = [ [100, 8, 3, 2, 3, 49], [38, 2, 20] ] # Mock fluorescence data
102
  return all_paths, path_lengths, measured_trace_fluorescence
103
 
104
  def test_process_cell_traces_return_type(mock_data):
 
121
  def test_process_cell_traces_peaks(mock_data):
122
  all_paths, path_lengths, measured_trace_fluorescence = mock_data
123
  result = process_cell_traces(all_paths, path_lengths, measured_trace_fluorescence)
124
+ print(result)
125
  peaks = [p.peaks for p in result.pathdata_list]
126
  assert peaks == [[0,5],[]]
127
 
128
  # Mock data
129
  @pytest.fixture
130
  def mock_celldata():
131
+ pathdata1 = PathData(peaks=[0, 5], points=[(0,0,0), (0,2,0), (0,5,0), (0,10,0), (0,15,0), (0,20,0)], removed_peaks=[], o_hei10=[100, 8, 3, 2, 3, 69], SC_length=2.2)
132
+ pathdata2 = PathData(peaks=[2], points=[(1,20,0), (1,20,10), (1,20,20) ], removed_peaks=[RemovedPeakData(0, (0,5))], o_hei10=[38, 2, 20], SC_length=2.3)
133
  return CellData(pathdata_list=[pathdata1, pathdata2])
134
 
135
+ def test_analyse_celldata(mock_celldata):
136
+ data_frame, foci_absolute_intensity, foci_position_index, dominated_foci_data, trace_median_intensity, trace_thresholds = analyse_celldata(mock_celldata, {'peak_threshold': 0.4, 'threshold_type':'per-trace'})
137
+ assert len(data_frame) == len(mock_celldata.pathdata_list), "Mismatch in dataframe length"
138
+ assert len(foci_absolute_intensity) == len(mock_celldata.pathdata_list), "Mismatch in relative intensities length"
139
+ assert len(foci_position_index) == len(mock_celldata.pathdata_list), "Mismatch in positions length"
 
140
 
141
+ assert list(map(list, foci_position_index)) == [[0, 5], [2]]
142
 
143
 
144
+ def test_analyse_celldata_per_cell(mock_celldata):
145
+ data_frame, foci_absolute_intensity, foci_position_index, dominated_foci_data, trace_median_intensity, trace_thresholds = analyse_celldata(mock_celldata, {'peak_threshold': 0.4, 'threshold_type':'per-cell'})
146
+ assert len(data_frame) == len(mock_celldata.pathdata_list), "Mismatch in relative intensities length"
147
+ assert len(foci_absolute_intensity) == len(mock_celldata.pathdata_list), "Mismatch in positions length"
148
+ assert len(foci_position_index) == len(mock_celldata.pathdata_list), "Mismatch in position indices length"
149
+ assert list(map(list, foci_position_index)) == [[0, 5], []]
150