Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Actin compression plots visualization #26

Merged
merged 6 commits into from
Dec 21, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
106 changes: 9 additions & 97 deletions examples/actin/docker/src/actin.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,23 +3,17 @@

import os
import argparse
import time

import numpy as np
import pandas
import psutil

from subcell_analysis.readdy import (
ReaddyLoader,
ReaddyPostProcessor,
)
from simulariumio import BinaryWriter

from simularium_readdy_models.actin import (
FiberData,
ActinSimulation,
ActinGenerator,
ActinTestData,
ActinStructure,
)
from simularium_readdy_models.visualization import ActinVisualization
from simularium_readdy_models import ReaddyUtil
Expand Down Expand Up @@ -114,98 +108,10 @@ def report_hardware_usage():
)


def analyze_results(parameters, save_pickle=False):
# get analysis parameters
plot_actin_compression = parameters.get("plot_actin_compression", False)
visualize_edges = parameters.get("visualize_edges", False)
visualize_normals = parameters.get("visualize_normals", False)
visualize_control_pts = parameters.get("visualize_control_pts", False)

# convert to simularium
traj_data = ActinVisualization.simularium_trajectory(
path_to_readdy_h5=parameters["name"] + ".h5",
box_size=parameters["box_size"],
total_steps=parameters["total_steps"],
time_multiplier=1e-3, # assume 1e3 recorded steps
longitudinal_bonds=bool(parameters.get("longitudinal_bonds", True)),
)

# load different views of ReaDDy data
post_processor = None
fiber_chain_ids = None
axis_positions = None
new_chain_ids = None
if visualize_normals or visualize_control_pts or plot_actin_compression or visualize_edges:
periodic_boundary = parameters.get("periodic_boundary", False)
post_processor = ReaddyPostProcessor(
trajectory=ReaddyLoader(
h5_file_path=parameters["name"] + ".h5",
min_time_ix=0,
max_time_ix=-1,
time_inc=1,
timestep=100.0,
save_pickle_file=save_pickle,
).trajectory(),
box_size=parameters["box_size"],
periodic_boundary=periodic_boundary,
)
if visualize_normals or visualize_control_pts or plot_actin_compression:
fiber_chain_ids = post_processor.linear_fiber_chain_ids(
start_particle_phrases=["pointed"],
other_particle_types=[
"actin#",
"actin#ATP_",
"actin#mid_",
"actin#mid_ATP_",
"actin#fixed_",
"actin#fixed_ATP_",
"actin#mid_fixed_",
"actin#mid_fixed_ATP_",
"actin#barbed_",
"actin#barbed_ATP_",
"actin#fixed_barbed_",
"actin#fixed_barbed_ATP_",
],
polymer_number_range=5,
)
if visualize_normals or visualize_control_pts:
axis_positions, new_chain_ids = post_processor.linear_fiber_axis_positions(
fiber_chain_ids=fiber_chain_ids,
ideal_positions=ActinStructure.mother_positions[2:5],
ideal_vector_to_axis=ActinStructure.vector_to_axis(),
)

# create plots
if plot_actin_compression:
print("plot actin compression")
traj_data.plots = ActinVisualization.generate_actin_compression_plots(
post_processor,
fiber_chain_ids,
temperature_c=parameters["temperature_C"],
)

# add annotation objects to the spatial data
traj_data = ActinVisualization.add_spatial_annotations(
traj_data,
post_processor,
visualize_edges,
visualize_normals,
visualize_control_pts,
new_chain_ids,
axis_positions,
)

# save simularium file
BinaryWriter.save(
trajectory_data=traj_data,
output_path=parameters["name"] + ".h5",
validate_ids=False,
)


def main():
args = parse_args()
parameters = setup_parameters(args)
start_time = time.time()
actin_simulation = ActinSimulation(
parameters=parameters,
record=True,
Expand All @@ -217,8 +123,14 @@ def main():
timestep=actin_simulation.parameters.get("internal_timestep", 0.1),
show_summary=False,
)
print("Run time: %s seconds " % (time.time() - start_time))
report_hardware_usage()
analyze_results(actin_simulation.parameters, args.save_pickle)
ActinVisualization.analyze_and_visualize_trajectory(
actin_simulation.parameters["name"],
actin_simulation.parameters["total_steps"],
actin_simulation.parameters,
args.save_pickle,
)


if __name__ == "__main__":
Expand Down
Binary file modified examples/actin/template.xlsx
Binary file not shown.
3 changes: 1 addition & 2 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -144,10 +144,9 @@ ignore = [
".cookiecutter.yaml",
"*docs/*",
"*.ipynb",
"*examples/*",
"*examples/**",
"*tests/**",
"environment.yml",
"examples/*",
]

[tool.mypy]
Expand Down
52 changes: 24 additions & 28 deletions simularium_readdy_models/actin/actin_analyzer.py
Original file line number Diff line number Diff line change
Expand Up @@ -1149,7 +1149,7 @@ def analyze_filament_length(
return np.array(filament_length)

@staticmethod
def analyze_bond_stretch(monomer_data, box_size, periodic_boundary, stride=1):
def analyze_bond_stretch(trajectory, box_size, periodic_boundary, stride=1):
"""
Get the difference in bond length from ideal
for lateral and longitudinal actin bonds.
Expand All @@ -1163,22 +1163,18 @@ def analyze_bond_stretch(monomer_data, box_size, periodic_boundary, stride=1):
ActinStructure.mother_positions[2] - ActinStructure.mother_positions[0]
)
print("Analyzing bond stretch...")
for time_index in range(0, len(monomer_data), stride):
for time_index in range(0, len(trajectory), stride):
stretch_lat.append([])
stretch_long.append([])
new_time_index = math.floor(time_index / stride)
filament = monomer_data[time_index]["topologies"][0]["particle_ids"]
filament = trajectory[time_index].topologies[0].particle_ids
for index in range(len(filament) - 2):
particle = monomer_data[time_index]["particles"][filament[index]]
particle_lat = monomer_data[time_index]["particles"][
filament[index + 1]
]
particle_long = monomer_data[time_index]["particles"][
filament[index + 2]
]
type_name = particle["type_name"]
type_name_lat = particle_lat["type_name"]
type_name_long = particle_long["type_name"]
particle = trajectory[time_index].particles[filament[index]]
particle_lat = trajectory[time_index].particles[filament[index + 1]]
particle_long = trajectory[time_index].particles[filament[index + 2]]
type_name = particle.type_name
type_name_lat = particle_lat.type_name
type_name_long = particle_long.type_name
if (
"fixed" in type_name
or "fixed" in type_name_lat
Expand All @@ -1187,9 +1183,9 @@ def analyze_bond_stretch(monomer_data, box_size, periodic_boundary, stride=1):
stretch_lat[new_time_index].append(0.0)
stretch_long[new_time_index].append(0.0)
continue
pos = particle["position"]
pos_lat = particle_lat["position"]
pos_long = particle_long["position"]
pos = particle.position
pos_lat = particle_lat.position
pos_long = particle_long.position
if periodic_boundary:
pos_lat = ReaddyUtil.get_non_periodic_boundary_position(
pos, pos_lat, box_size
Expand All @@ -1208,7 +1204,7 @@ def analyze_bond_stretch(monomer_data, box_size, periodic_boundary, stride=1):
return np.array(stretch_lat), np.array(stretch_long)

@staticmethod
def analyze_angle_stretch(monomer_data, box_size, periodic_boundary, stride=1):
def analyze_angle_stretch(trajectory, box_size, periodic_boundary, stride=1):
"""
Get the difference in angles (degrees) from ideal
for actin angles between:
Expand All @@ -1223,24 +1219,24 @@ def analyze_angle_stretch(monomer_data, box_size, periodic_boundary, stride=1):
ideal_angle_lat_long = ActinStructure.actin_to_actin_angle(True, False, True)
ideal_angle_long_long = ActinStructure.actin_to_actin_angle(False, False, True)
print("Analyzing angle stretch...")
for time_index in range(0, len(monomer_data), stride):
for time_index in range(0, len(trajectory), stride):
result_lat_lat.append([])
result_lat_long.append([])
result_long_long.append([])
new_time_index = math.floor(time_index / stride)
filament = monomer_data[time_index]["topologies"][0]["particle_ids"]
filament = trajectory[time_index].topologies[0].particle_ids
for index in range(len(filament) - 4):
particles = []
positions = []
fixed = False
for d in range(5):
particles.append(
monomer_data[time_index]["particles"][filament[index + d]]
trajectory[time_index].particles[filament[index + d]]
)
if "fixed" in particles[d]["type_name"]:
if "fixed" in particles[d].type_name:
fixed = True
break
positions.append(particles[d]["position"])
positions.append(particles[d].position)
if d > 0 and periodic_boundary:
positions[d] = ReaddyUtil.get_non_periodic_boundary_position(
positions[d - 1], positions[d], box_size
Expand Down Expand Up @@ -1287,7 +1283,7 @@ def analyze_angle_stretch(monomer_data, box_size, periodic_boundary, stride=1):
)

@staticmethod
def analyze_dihedral_stretch(monomer_data, box_size, periodic_boundary, stride=1):
def analyze_dihedral_stretch(trajectory, box_size, periodic_boundary, stride=1):
"""
Get the difference in dihedral angles (degrees) from ideal
for actin angles between:
Expand All @@ -1305,23 +1301,23 @@ def analyze_dihedral_stretch(monomer_data, box_size, periodic_boundary, stride=1
False, False, False, True
)
print("Analyzing dihedral stretch...")
for time_index in range(0, len(monomer_data), stride):
for time_index in range(0, len(trajectory), stride):
result_lat_lat_lat.append([])
result_long_long_long.append([])
new_time_index = math.floor(time_index / stride)
filament = monomer_data[time_index]["topologies"][0]["particle_ids"]
filament = trajectory[time_index].topologies[0].particle_ids
for index in range(len(filament) - 6):
particles = []
positions = []
fixed = False
for d in range(7):
particles.append(
monomer_data[time_index]["particles"][filament[index + d]]
trajectory[time_index].particles[filament[index + d]]
)
if "fixed" in particles[d]["type_name"]:
if "fixed" in particles[d].type_name:
fixed = True
break
positions.append(particles[d]["position"])
positions.append(particles[d].position)
if d > 0 and periodic_boundary:
positions[d] = ReaddyUtil.get_non_periodic_boundary_position(
positions[d - 1], positions[d], box_size
Expand Down
Loading