Coverage for subcell_pipeline/analysis/compression_metrics/polymer_trace.py: 17%
146 statements
« prev ^ index » next coverage.py v7.5.3, created at 2024-08-29 15:14 +0000
« prev ^ index » next coverage.py v7.5.3, created at 2024-08-29 15:14 +0000
1"""Methods to calculate metrics from polymer trace data."""
3from typing import Any
5import numpy as np
6from sklearn.decomposition import PCA
8from subcell_pipeline.analysis.compression_metrics.constants import (
9 ABSOLUTE_TOLERANCE,
10 DEFAULT_BENDING_CONSTANT,
11)
12from subcell_pipeline.analysis.compression_metrics.vectors import (
13 get_end_to_end_unit_vector,
14)
17def get_end_to_end_axis_distances_and_projections(
18 polymer_trace: np.ndarray,
19) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
20 """
21 Calculate distances of the polymer trace points from the end-to-end axis.
23 Here, the end-to-end axis is defined as the line joining the first and last
24 fiber point.
26 Parameters
27 ----------
28 polymer_trace
29 Array containing the x,y,z positions of the polymer trace.
31 Returns
32 -------
33 perp_distances
34 Perpendicular distances of the polymer trace from the end-to-end axis.
35 scaled_projections
36 Length of fiber point projections along the end-to-end axis, scaled by
37 axis length. Can be negative.
38 projection_positions
39 Positions of points on the end-to-end axis that are closest from the
40 respective points in the polymer trace. The distance from
41 projection_positions to the trace points is the shortest distance from
42 the end-to-end axis.
43 """
44 end_to_end_axis = get_end_to_end_unit_vector(polymer_trace=polymer_trace)
45 end_to_end_axis_length = np.linalg.norm(polymer_trace[-1] - polymer_trace[0])
47 position_vectors = polymer_trace - polymer_trace[0]
48 projections = np.dot(position_vectors, end_to_end_axis)
49 scaled_projections = projections / end_to_end_axis_length
51 projection_positions = polymer_trace[0] + projections[:, None] * end_to_end_axis
53 perp_distances = np.linalg.norm(polymer_trace - projection_positions, axis=1)
55 return perp_distances, scaled_projections, projection_positions
58def get_average_distance_from_end_to_end_axis(
59 polymer_trace: np.ndarray, **options: dict[str, Any]
60) -> float:
61 """
62 Calculate the average perpendicular distance of polymer trace points from
63 the end-to-end axis.
65 Parameters
66 ----------
67 polymer_trace
68 Array containing the x,y,z positions of the polymer trace.
69 **options
70 Additional options as key-value pairs. Unused.
72 Returns
73 -------
74 :
75 Average perpendicular distance of polymer trace points from the
76 end-to-end axis.
77 """
78 perp_distances, _, _ = get_end_to_end_axis_distances_and_projections(
79 polymer_trace=polymer_trace
80 )
81 avg_perp_distance = np.nanmean(perp_distances)
83 return float(avg_perp_distance)
86def get_asymmetry_of_peak(
87 polymer_trace: np.ndarray, **options: dict[str, Any]
88) -> float:
89 """
90 Calculate the scaled distance of the projection of the peak from the
91 end-to-end axis midpoint.
93 Parameters
94 ----------
95 polymer_trace
96 Array containing the x,y,z positions of the polymer trace.
97 **options
98 Additional options as key-value pairs. Unused.
100 Returns
101 -------
102 :
103 Scaled distance of the projection of the peak from the axis midpoint.
104 """
105 (
106 perp_distances,
107 scaled_projections,
108 _,
109 ) = get_end_to_end_axis_distances_and_projections(polymer_trace=polymer_trace)
111 # if all perpendicular distances are zero, return 0
112 if np.all(perp_distances < ABSOLUTE_TOLERANCE):
113 return 0
115 projection_of_peak = scaled_projections[perp_distances == np.max(perp_distances)]
116 peak_asym = np.max(
117 np.abs(projection_of_peak - 0.5)
118 ) # max kinda handles multiple peaks
120 return peak_asym
123def get_pca_polymer_trace_projection(polymer_trace: np.ndarray) -> np.ndarray:
124 """
125 Calculate the PCA projection of the polymer trace.
127 Parameters
128 ----------
129 polymer_trace
130 Array containing the x,y,z positions of the polymer trace.
132 Returns
133 -------
134 :
135 PCA projection of the polymer trace.
136 """
137 pca = fit_pca_to_polymer_trace(polymer_trace=polymer_trace)
138 return pca.transform(polymer_trace)
141def get_contour_length_from_trace(
142 polymer_trace: np.ndarray, **options: dict[str, Any]
143) -> float:
144 """
145 Calculate the sum of inter-monomer distances in the trace.
147 Parameters
148 ----------
149 polymer_trace
150 Array containing the x,y,z positions of the polymer trace.
151 **options
152 Additional options as key-value pairs. Unused.
154 Returns
155 -------
156 :
157 Sum of inter-monomer distances in the trace.
158 """
159 return np.sum(np.linalg.norm(np.diff(polymer_trace, axis=0), axis=1))
162def get_bending_energy_from_trace(
163 polymer_trace: np.ndarray, **options: dict[str, Any]
164) -> float:
165 """
166 Calculate the bending energy per monomer of a polymer trace.
168 Parameters
169 ----------
170 polymer_trace
171 Array containing the x,y,z positions of the polymer trace.
172 **options
173 Additional options as key-value pairs:
175 bending_constant: float
176 Bending constant of the fiber in pN nm.
178 Returns
179 -------
180 :
181 Bending energy per monomer of the polymer trace.
182 """
183 bending_constant = options.get("bending_constant", DEFAULT_BENDING_CONSTANT)
185 assert isinstance(bending_constant, (float, np.floating))
187 cos_angle = np.zeros(len(polymer_trace) - 2)
188 for ind in range(len(polymer_trace) - 2):
189 vec1 = polymer_trace[ind + 1] - polymer_trace[ind]
190 vec2 = polymer_trace[ind + 2] - polymer_trace[ind + 1]
192 if np.isclose(np.linalg.norm(vec1), 0.0) or np.isclose(
193 np.linalg.norm(vec2), 0.0
194 ):
195 # TODO handle this differently?
196 cos_angle[ind] = 0.0
197 print("Warning: zero vector in bending energy calculation.")
198 continue
200 cos_angle[ind] = (
201 np.dot(vec1, vec2) / np.linalg.norm(vec1) / np.linalg.norm(vec2)
202 )
204 # since the bending constant is obtained from a kwargs dictionary
205 # the type checker is unable to infer its type
206 energy = bending_constant * (1 - np.nanmean(cos_angle))
208 return energy.item()
211def get_2d_polymer_trace(
212 polymer_trace: np.ndarray, compression_axis: int = 0
213) -> np.ndarray:
214 """
215 Get the 2D projection of the polymer trace.
217 Parameters
218 ----------
219 polymer_trace
220 Array containing the x,y,z positions of the polymer trace.
221 compression_axis
222 The axis along which the polymer trace is compressed.
224 Returns
225 -------
226 :
227 The 2D projection of the polymer trace.
228 """
229 if polymer_trace.shape[1] == 2:
230 return polymer_trace
232 return polymer_trace[
233 :, [ax for ax in range(polymer_trace.shape[1]) if ax != compression_axis]
234 ]
237def get_normalized_tangent_vectors(polymer_trace: np.ndarray) -> np.ndarray:
238 """
239 Calculate the normalized tangent vectors for a polymer trace.
241 Parameters
242 ----------
243 polymer_trace
244 Array containing the x,y,z positions of the polymer trace.
246 Returns
247 -------
248 :
249 The normalized tangent vectors for the polymer.
250 """
251 tangents = np.diff(polymer_trace, axis=0)
253 tangent_lengths = np.linalg.norm(tangents, axis=1)
255 if np.all(tangent_lengths < ABSOLUTE_TOLERANCE):
256 return np.zeros_like(tangents)
258 tangents /= tangent_lengths[:, np.newaxis]
260 return tangents
263def get_twist_angle(polymer_trace: np.ndarray, **options: dict[str, Any]) -> float:
264 """
265 Calculate the twist angle of the polymer trace.
267 Parameters
268 ----------
269 polymer_trace
270 Array containing the x,y,z positions of the polymer trace.
271 **options
272 Additional options as key-value pairs. Unused.
274 Returns
275 -------
276 :
277 Twist angle of the polymer trace in degrees.
278 """
279 compression_axis = options.get("compression_axis", 0)
280 assert isinstance(compression_axis, int)
282 trace_2d = get_2d_polymer_trace(
283 polymer_trace=polymer_trace, compression_axis=compression_axis
284 )
286 tangents = get_normalized_tangent_vectors(
287 polymer_trace=trace_2d,
288 )
290 angle = get_angle_between_vectors(tangents[0], -tangents[-1], signed=False)
291 chirality = get_chirality(polymer_trace=polymer_trace)
293 return chirality * angle * 180 / np.pi
296def get_chirality(polymer_trace: np.ndarray, **options: dict[str, Any]) -> float:
297 """
298 Calculate the chirality of a polymer trace.
300 Parameters
301 ----------
302 polymer_trace
303 Array containing the x,y,z positions of the polymer trace.
304 **options
305 Additional options as key-value pairs. Unused.
307 Returns
308 -------
309 :
310 Chirality of the polymer trace.
311 """
312 trace_2d = get_2d_polymer_trace(polymer_trace=polymer_trace)
313 tangents = get_normalized_tangent_vectors(polymer_trace=trace_2d)
315 chirality = 0
316 for i in range(len(tangents) - 1):
317 cross_product = np.cross(tangents[i], tangents[i + 1])
318 if cross_product > 0:
319 chirality += 1
320 elif cross_product < 0:
321 chirality -= 1
323 return np.sign(chirality)
326def get_total_fiber_twist(
327 polymer_trace: np.ndarray, **options: dict[str, Any]
328) -> float:
329 """
330 Calculate the total twist of a polymer trace using the normal vectors.
332 Parameters
333 ----------
334 polymer_trace
335 Array containing the x,y,z positions of the polymer trace.
336 **options
337 Additional options as key-value pairs:
339 signed: bool
340 Whether to return the signed or unsigned total twist.
342 Returns
343 -------
344 :
345 Total twist of the polymer trace.
346 """
347 signed = options.get("signed", False)
348 assert isinstance(signed, bool)
350 tangents = np.diff(polymer_trace, axis=0)
351 tangents /= np.linalg.norm(tangents, axis=1)[:, np.newaxis]
353 normals = np.cross(tangents[:-1], tangents[1:])
354 normal_lengths = np.linalg.norm(normals, axis=1)
355 if np.all(normal_lengths < ABSOLUTE_TOLERANCE):
356 return 0
357 normals = normals / normal_lengths[:, np.newaxis]
359 twists = []
360 for i in range(1, len(normals)):
361 angle = get_angle_between_vectors(normals[i - 1], normals[i], signed=signed)
363 twists.append(angle)
365 # Sum the twist angles to get the total twist
366 total_twist = np.sum(twists)
368 # Normalize by the contour length
369 total_twist /= get_contour_length_from_trace(polymer_trace)
371 return total_twist
374def get_total_fiber_twist_project(
375 polymer_trace: np.ndarray, **options: dict[str, Any]
376) -> float:
377 """
378 Calculate the total twist using projections of the polymer trace in the 2nd
379 and 3rd dimension.
381 Parameters
382 ----------
383 polymer_trace
384 Array containing the x,y,z positions of the polymer trace.
385 **options
386 Additional options as key-value pairs:
388 compression_axis: int
389 Axis along which the polymer trace is compressed.
390 signed: bool
391 Whether to return the signed or unsigned total twist.
392 tolerance: float
393 Tolerance for vector length.
395 Returns
396 -------
397 :
398 Sum of angles between PCA projection vectors.
399 """
400 compression_axis = options.get("compression_axis", 0)
401 signed = options.get("signed", True)
402 tolerance = options.get("tolerance", ABSOLUTE_TOLERANCE)
404 assert isinstance(signed, bool)
405 assert isinstance(tolerance, (float, np.floating))
406 assert isinstance(compression_axis, int)
408 trace_2d = get_2d_polymer_trace(
409 polymer_trace=polymer_trace, compression_axis=compression_axis
410 )
411 trace_2d = trace_2d - np.mean(trace_2d, axis=0)
413 return get_total_fiber_twist_2d(trace_2d, signed=signed, tolerance=tolerance)
416def get_total_fiber_twist_pca(
417 polymer_trace: np.ndarray, tolerance: float = ABSOLUTE_TOLERANCE
418) -> float:
419 """
420 Calculate the total twist using PCA projections of the polymer trace in the
421 2nd and 3rd dimension.
423 Parameters
424 ----------
425 polymer_trace
426 Array containing the x,y,z positions of the polymer trace.
427 tolerance
428 Tolerance for vector length.
430 Returns
431 -------
432 :
433 Sum of angles between PCA projection vectors.
434 """
435 pca_trace = get_pca_polymer_trace_projection(polymer_trace=polymer_trace)
436 pca_trace_2d = pca_trace[:, 1:]
438 return get_total_fiber_twist_2d(pca_trace_2d, tolerance=tolerance)
441def get_angle_between_vectors(
442 vec1: np.ndarray, vec2: np.ndarray, signed: bool = False
443) -> float:
444 """
445 Calculate the signed angle between two vectors.
447 Parameters
448 ----------
449 vec1
450 The first vector.
451 vec2
452 The second vector.
453 signed
454 True to get signed angle between vectors, False otherwise.
456 Returns
457 -------
458 :
459 Angle between vec1 and vec2.
460 """
461 vec1_length = np.linalg.norm(vec1)
462 vec2_length = np.linalg.norm(vec2)
464 if vec1_length < ABSOLUTE_TOLERANCE or vec2_length < ABSOLUTE_TOLERANCE:
465 return 0
467 vec1 = vec1 / vec1_length
468 vec2 = vec2 / vec2_length
470 angle = np.arccos(np.dot(vec1, vec2))
472 if signed:
473 if np.cross(vec1, vec2) < 0:
474 angle = -angle
476 return angle
479def get_total_fiber_twist_2d(
480 trace_2d: np.ndarray, signed: bool = False, tolerance: float = ABSOLUTE_TOLERANCE
481) -> float:
482 """
483 Calculate the total twist for 2D traces.
485 The 2D twist is defined as the sum of (signed) angles between consecutive
486 vectors in the 2D projection along the compression axis.
488 Parameters
489 ----------
490 trace_2d
491 Array containing the x,y positions of the polymer trace.
492 signed
493 True to calculate signed total twist, False otherwise.
494 tolerance
495 Tolerance for vector length.
497 Returns
498 -------
499 :
500 Sum of angles between trace vectors.
501 """
502 prev_vec = None
503 angles = np.zeros(len(trace_2d))
504 for i in range(len(trace_2d)):
505 if prev_vec is None:
506 prev_vec_length = np.linalg.norm(trace_2d[i])
507 if prev_vec_length < tolerance:
508 prev_vec = None
509 continue
510 prev_vec = trace_2d[i] / prev_vec_length
512 curr_vec_length = np.linalg.norm(trace_2d[i])
513 if curr_vec_length < tolerance:
514 continue
515 curr_vec = trace_2d[i] / curr_vec_length
517 angles[i] = get_angle_between_vectors(prev_vec, curr_vec, signed=signed)
519 prev_vec = curr_vec
521 return np.abs(np.nansum(angles) / 2 / np.pi)
524def fit_pca_to_polymer_trace(polymer_trace: np.ndarray) -> PCA:
525 """
526 Fits PCA to the polymer trace and returns the fitted object.
528 Parameters
529 ----------
530 polymer_trace
531 Array containing the x,y,z positions of the polymer trace.
533 Returns
534 -------
535 :
536 PCA object fit to the polymer trace.
537 """
538 pca = PCA(n_components=3)
539 pca.fit(polymer_trace)
540 return pca
543def get_third_component_variance(
544 polymer_trace: np.ndarray, **options: dict[str, Any]
545) -> float:
546 """
547 Calculate the 3rd PCA component given the x,y,z positions of a fiber.
549 This component reflects non-coplanarity/out-of-planeness of the fiber.
551 Parameters
552 ----------
553 polymer_trace
554 Array containing the x,y,z positions of the polymer trace.
555 **options
556 Additional options as key-value pairs. Unused.
558 Returns
559 -------
560 :
561 Variance explained by the third principal component.
562 """
563 pca = fit_pca_to_polymer_trace(polymer_trace=polymer_trace)
564 return pca.explained_variance_ratio_[2]
567def get_sum_bending_energy(
568 fiber_energy: np.ndarray, **options: dict[str, Any]
569) -> float:
570 """
571 Calculate the sum of bending energy from the given fiber energy array.
573 Parameters
574 ----------
575 fiber_energy
576 Array containing fiber energy values.
577 **options
578 Additional options as key-value pairs. Unused.
580 Returns
581 -------
582 :
583 Sum of bending energy.
584 """
585 return fiber_energy[3].sum()
588def get_compression_ratio(
589 polymer_trace: np.ndarray, **options: dict[str, Any]
590) -> float:
591 """
592 Calculate the compression ratio of a polymer trace.
594 The compression ratio is defined as 1 minus the ratio of the length of the
595 end-to-end vector to the contour length of the polymer trace.
597 Parameters
598 ----------
599 polymer_trace
600 Array containing the x,y,z positions of the polymer trace.
601 **options
602 Additional options as key-value pairs. Unused.
604 Returns
605 -------
606 :
607 The compression ratio of the polymer trace.
608 """
609 end_to_end_axis_length = np.linalg.norm(polymer_trace[-1] - polymer_trace[0]).item()
610 return 1 - end_to_end_axis_length / get_contour_length_from_trace(polymer_trace)