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

1"""Methods to calculate metrics from polymer trace data.""" 

2 

3from typing import Any 

4 

5import numpy as np 

6from sklearn.decomposition import PCA 

7 

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) 

15 

16 

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. 

22 

23 Here, the end-to-end axis is defined as the line joining the first and last 

24 fiber point. 

25 

26 Parameters 

27 ---------- 

28 polymer_trace 

29 Array containing the x,y,z positions of the polymer trace. 

30 

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]) 

46 

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 

50 

51 projection_positions = polymer_trace[0] + projections[:, None] * end_to_end_axis 

52 

53 perp_distances = np.linalg.norm(polymer_trace - projection_positions, axis=1) 

54 

55 return perp_distances, scaled_projections, projection_positions 

56 

57 

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. 

64 

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. 

71 

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) 

82 

83 return float(avg_perp_distance) 

84 

85 

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. 

92 

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. 

99 

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) 

110 

111 # if all perpendicular distances are zero, return 0 

112 if np.all(perp_distances < ABSOLUTE_TOLERANCE): 

113 return 0 

114 

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 

119 

120 return peak_asym 

121 

122 

123def get_pca_polymer_trace_projection(polymer_trace: np.ndarray) -> np.ndarray: 

124 """ 

125 Calculate the PCA projection of the polymer trace. 

126 

127 Parameters 

128 ---------- 

129 polymer_trace 

130 Array containing the x,y,z positions of the polymer trace. 

131 

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) 

139 

140 

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. 

146 

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. 

153 

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)) 

160 

161 

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. 

167 

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: 

174 

175 bending_constant: float 

176 Bending constant of the fiber in pN nm. 

177 

178 Returns 

179 ------- 

180 : 

181 Bending energy per monomer of the polymer trace. 

182 """ 

183 bending_constant = options.get("bending_constant", DEFAULT_BENDING_CONSTANT) 

184 

185 assert isinstance(bending_constant, (float, np.floating)) 

186 

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] 

191 

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 

199 

200 cos_angle[ind] = ( 

201 np.dot(vec1, vec2) / np.linalg.norm(vec1) / np.linalg.norm(vec2) 

202 ) 

203 

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)) 

207 

208 return energy.item() 

209 

210 

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. 

216 

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. 

223 

224 Returns 

225 ------- 

226 : 

227 The 2D projection of the polymer trace. 

228 """ 

229 if polymer_trace.shape[1] == 2: 

230 return polymer_trace 

231 

232 return polymer_trace[ 

233 :, [ax for ax in range(polymer_trace.shape[1]) if ax != compression_axis] 

234 ] 

235 

236 

237def get_normalized_tangent_vectors(polymer_trace: np.ndarray) -> np.ndarray: 

238 """ 

239 Calculate the normalized tangent vectors for a polymer trace. 

240 

241 Parameters 

242 ---------- 

243 polymer_trace 

244 Array containing the x,y,z positions of the polymer trace. 

245 

246 Returns 

247 ------- 

248 : 

249 The normalized tangent vectors for the polymer. 

250 """ 

251 tangents = np.diff(polymer_trace, axis=0) 

252 

253 tangent_lengths = np.linalg.norm(tangents, axis=1) 

254 

255 if np.all(tangent_lengths < ABSOLUTE_TOLERANCE): 

256 return np.zeros_like(tangents) 

257 

258 tangents /= tangent_lengths[:, np.newaxis] 

259 

260 return tangents 

261 

262 

263def get_twist_angle(polymer_trace: np.ndarray, **options: dict[str, Any]) -> float: 

264 """ 

265 Calculate the twist angle of the polymer trace. 

266 

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. 

273 

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) 

281 

282 trace_2d = get_2d_polymer_trace( 

283 polymer_trace=polymer_trace, compression_axis=compression_axis 

284 ) 

285 

286 tangents = get_normalized_tangent_vectors( 

287 polymer_trace=trace_2d, 

288 ) 

289 

290 angle = get_angle_between_vectors(tangents[0], -tangents[-1], signed=False) 

291 chirality = get_chirality(polymer_trace=polymer_trace) 

292 

293 return chirality * angle * 180 / np.pi 

294 

295 

296def get_chirality(polymer_trace: np.ndarray, **options: dict[str, Any]) -> float: 

297 """ 

298 Calculate the chirality of a polymer trace. 

299 

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. 

306 

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) 

314 

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 

322 

323 return np.sign(chirality) 

324 

325 

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. 

331 

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: 

338 

339 signed: bool 

340 Whether to return the signed or unsigned total twist. 

341 

342 Returns 

343 ------- 

344 : 

345 Total twist of the polymer trace. 

346 """ 

347 signed = options.get("signed", False) 

348 assert isinstance(signed, bool) 

349 

350 tangents = np.diff(polymer_trace, axis=0) 

351 tangents /= np.linalg.norm(tangents, axis=1)[:, np.newaxis] 

352 

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] 

358 

359 twists = [] 

360 for i in range(1, len(normals)): 

361 angle = get_angle_between_vectors(normals[i - 1], normals[i], signed=signed) 

362 

363 twists.append(angle) 

364 

365 # Sum the twist angles to get the total twist 

366 total_twist = np.sum(twists) 

367 

368 # Normalize by the contour length 

369 total_twist /= get_contour_length_from_trace(polymer_trace) 

370 

371 return total_twist 

372 

373 

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. 

380 

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: 

387 

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. 

394 

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) 

403 

404 assert isinstance(signed, bool) 

405 assert isinstance(tolerance, (float, np.floating)) 

406 assert isinstance(compression_axis, int) 

407 

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) 

412 

413 return get_total_fiber_twist_2d(trace_2d, signed=signed, tolerance=tolerance) 

414 

415 

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. 

422 

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. 

429 

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:] 

437 

438 return get_total_fiber_twist_2d(pca_trace_2d, tolerance=tolerance) 

439 

440 

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. 

446 

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. 

455 

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) 

463 

464 if vec1_length < ABSOLUTE_TOLERANCE or vec2_length < ABSOLUTE_TOLERANCE: 

465 return 0 

466 

467 vec1 = vec1 / vec1_length 

468 vec2 = vec2 / vec2_length 

469 

470 angle = np.arccos(np.dot(vec1, vec2)) 

471 

472 if signed: 

473 if np.cross(vec1, vec2) < 0: 

474 angle = -angle 

475 

476 return angle 

477 

478 

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. 

484 

485 The 2D twist is defined as the sum of (signed) angles between consecutive 

486 vectors in the 2D projection along the compression axis. 

487 

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. 

496 

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 

511 

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 

516 

517 angles[i] = get_angle_between_vectors(prev_vec, curr_vec, signed=signed) 

518 

519 prev_vec = curr_vec 

520 

521 return np.abs(np.nansum(angles) / 2 / np.pi) 

522 

523 

524def fit_pca_to_polymer_trace(polymer_trace: np.ndarray) -> PCA: 

525 """ 

526 Fits PCA to the polymer trace and returns the fitted object. 

527 

528 Parameters 

529 ---------- 

530 polymer_trace 

531 Array containing the x,y,z positions of the polymer trace. 

532 

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 

541 

542 

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. 

548 

549 This component reflects non-coplanarity/out-of-planeness of the fiber. 

550 

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. 

557 

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] 

565 

566 

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. 

572 

573 Parameters 

574 ---------- 

575 fiber_energy 

576 Array containing fiber energy values. 

577 **options 

578 Additional options as key-value pairs. Unused. 

579 

580 Returns 

581 ------- 

582 : 

583 Sum of bending energy. 

584 """ 

585 return fiber_energy[3].sum() 

586 

587 

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. 

593 

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. 

596 

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. 

603 

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)