Calculate arbitrary recentering operation from rotation matrix?

Cryo-EM Data Processing
1

Hi,

Consider a situation where we have non-point group symmetry in a molecule (combination of arbitrary rotation and translation). I would like to recenter based on alignment of two chains in a dimer - let’s say A & B of 7PTY in screenshot.

I can get the rotation matrix &/or axis/angle notation from Chimera, by using measure rotation after aligning with matchmaker - how can I convert this to ZYZ Euler angle notation with associated new center to use in volume alignment tools? Goal is to perform local refinement of both protomers in a dimer which does not have point group symmetry (but where the protomer itself is rigid enough to make local refinement useful).

Cheers
Oli

image
image1920×1092 201 KB

EDIT: The other potential issue here is that I would need to make a copy of the particle set with different particle UIDs - as is done in symmetry expansion, but without performing symmetry expansion

1 Like
2

Hi @olibclarke! This is an interesting question!

First off, I’ve recorded this as a feature request!

I’ll also take a look into a cs-tools implementation of this. It’s proving more challenging than I expected, but I’ll keep you posted if I figure it out.

In case you or others would like to take a swing at it, my general idea is to

  1. Load in the particles from a C1 symmetry expansion of your starting job. The C1 symmetry expansion doesn’t do anything except create the necessary fields we can use to mark the expanded particles as expanded copies.
  2. Programmatically create two Volume Align Tools jobs, one to rotate and one to shift. ChimeraX shifts after rotation, while CryoSPARC shifts first, so we need to do them in two steps.
  3. Finally, load in these particles and change their sym_expand/idx to 1, then save those as an external result.

The rotation matrix makes enough sense, but the translation part is proving more challenging. The ChimeraX docs say this about the rotation/translation matrix:

the first three columns of the matrix describe a rotation and the fourth describes a translation (performed after the rotation). The transformation is also described as an axis of rotation (a unit vector), point on the axis, degrees of rotation, and shift parallel to the axis

So, if we write a function to parse out these two components:

def parse_chimerax_matrix_string(chimx_mat:str) -> np.array:
    full_matrix = np.array([[float(y) for y in x.split(' ')] for x in chimx_mat.split('\n')])
    rot_matrix = full_matrix[:,:3]
    translation = full_matrix[:, 3].T
    return (rot_matrix, translation)

chimx_mat = """0.71327472 0.45534351 0.53282404 -94.10678621
-0.17014233 0.84996729 -0.49860525 111.25150489
-0.67991967 0.26498660 0.68373339 100.81847201"""
rot_mat, translation = parse_chimerax_matrix_string(chimx_mat)

we can take a look at each part. I’ll use the scipy library, which provides a nice wrapper for handling rotations (which I find very confusing to deal with on my own). To generate this example matrix, I rotated and aligned the T20S proteasome, so I should get an axis/angle rotation of around 51.4 degrees.

> from scipy.spatial.transform import Rotation as R
> rot_vec = R.from_matrix(rot_mat).as_rotvec()
> np.rad2deg(np.linalg.norm(rot_vec))
51.42872549572835

Looks good. We can then convert this to ZYZ Euler angles using the same library:

euler_angles = R.from_matrix(rot_mat).as_euler('zyz')

However, when I use these angles in a Volume Align Tools job, the results are clearly wrong:

image
image1214×858 91 KB

Maybe ChimeraX is rotating about the origin, which could also explain the massive shifts we have stored in translation: [296.87181427, 609.20380236, 593.33607243]. But we’ve basically reached the limit of my ability to abstractly model and understand 3D rotation in my head.

3 Likes
3

Ah, I mean to say as well: regarding the UIDs, you could solve this in two ways.

The easiest would be to combine it with the rotation step using the Reassign UIDs parameter in Volume Alignment Tools.

CryoSPARC tools can also re-assign UIDs with dataset.reassign_uids()

2 Likes
4

Thanks Rich this is very helpful! Chimera is not rotating around the origin I think, but instead around the calculated axis (indicated as a gray cylinder in the screenshot). I think this is also given as a vector in the log (along with a point with which the axis intersects)

1 Like
5

Indeed! But I couldn’t figure out why the rotation given by the matrix wasn’t producing the right orientation (regardless of any offset). After a bit away from the problem, I thought to double-check scipy’s documentation. “ZYZ” is intrinsic rotation (which we want); I had lazily typed “zyz” which is extrinsic :upside_down_face:.

The rest worked as expected. Script is below. A few tricky bits:

  • ChimeraX treats 0,0,0 as the origin, while CryoSPARC treats the middle of the box as the origin. So before performing the fitting, if you run volume #1 origin -C,-C,-C where C is box size / 2 * apix, the resulting rotation matrix will work right
  • If you want to be able to use sym_expand/idx to keep track of your virtual particles, make sure you’ve done a C1 symmetry expansion as I mentioned. If you don’t want to do that, you can skip the last bit and just use the job created by shifted_job.

You can directly copy-paste the ChimeraX matrix into the multi-line string, newlines and all.

from cryosparc.tools import CryoSPARC
import json
import numpy as np
from pathlib import Path
from scipy.spatial.transform import Rotation as R
with open(Path('~/instance-info.json').expanduser(), 'r') as f:
    instance_info = json.load(f)

cs = CryoSPARC(**instance_info)
assert cs.test_connection()

project_number = "P294"
workspace_number = "W2"

# src should be a symmetry expansion in C1 if you want
# to use the sym_expand/idx later
src_job_number = "J125"
src_job_particles_name = "particles"
src_job_vol_name = "volume"

lane_to_use = "cryoem10"

project = cs.find_project(project_number)
workspace = project.find_workspace(workspace_number)
job = project.find_job(src_job_number)
unrotated_particles = job.load_output(src_job_particles_name)

def parse_chimerax_matrix_string(chimx_mat:str) -> np.array:
    full_matrix = np.array(
        [
            [float(y) for y in x.split(' ')]
            for x in chimx_mat.split('\n')
        ]
    )
    rot_matrix = full_matrix[:,:3]
    translation = full_matrix[:, 3].T
    return (rot_matrix, translation)


chimx_mat = """0.71326635 -0.17010993 -0.67993657 -9.96570500
0.45530546 0.84999332 0.26496848 4.13432096
0.53286777 -0.49857193 0.68372361 -4.14917475"""
rot_mat, translation = parse_chimerax_matrix_string(chimx_mat)

# we want a new center in voxels, not a translation in A
translation /= unrotated_particles['alignments3D/psize_A'][0]
new_box_center = unrotated_particles['blob/shape'][0][0] / 2 - translation

euler = R.from_matrix(rot_mat).as_euler('ZYZ')

# chimeraX translates *after* rotating, CryoSPARC translates *before*.
# thus, we need two VAT jobs. One rotating, one translating.
unshifted_job = workspace.create_job(
    type = "volume_alignment_tools",
    connections = {
        "particles": (src_job_number, src_job_particles_name),
        "volume": (src_job_number, src_job_vol_name)
    },
    params = {
        # yes, this is the right parameter name for the Euler angles
        'recenter_axang': ','.join(str(x) for x in euler)
    }
)
unshifted_job.queue(lane_to_use)
unshifted_job.wait_for_done()

unshift_id = unshifted_job.doc['uid']
shifted_job = workspace.create_job(
    type = "volume_alignment_tools",
    connections = {
        "particles": (unshift_id, "particles"),
        "volume": (unshift_id, "volume")
    },
    params = {
        "recenter_shift": ",".join(str(x) for x in new_box_center)
    }
)
shifted_job.queue(lane_to_use)
shifted_job.wait_for_done()

shifted_id = shifted_job.doc['uid']
shifted_results = project.find_job(shifted_id)
shifted_particles = shifted_results.load_output('particles')
shifted_particles['sym_expand/idx'] = 1
project.save_external_result(
    workspace_uid=workspace_number,
    dataset=shifted_particles,
    type='particle',
    slots=['sym_expand'],
    passthrough=(shifted_id, 'particles'),
    title='non-point-group sym expand'
)

And now, for my usual caveat. I’ve tested this by

  1. arbitrarily rotating T20S proteasome to ensure that the symmetry axis is not aligned with Z
  2. manually performing one of the 7-fold rotations and using fitMap
  3. using the resulting matrix in this script

Here is the result: blue is the original map, orange the result of the final save_external_result:

image5
image51000×1000 87 KB

I think that should work for your use case, but let me know if not!

3 Likes
6

Awesome thanks, that should do the trick!!

1 Like
7

Hi Rich, I tried this, but it gave the following error:

(cryosparc-tools) user@ubuntu:~$ python realign_script.py
Connection succeeded to CryoSPARC command_core at http://localhost:39002
Connection succeeded to CryoSPARC command_vis at http://localhost:39003
Connection succeeded to CryoSPARC command_rtp at http://localhost:39005
Traceback (most recent call last):
  File "realign_script.py", line 43, in <module>
    rot_mat, translation = parse_chimerax_matrix_string(chimx_mat)
  File "realign_script.py", line 30, in parse_chimerax_matrix_string
    [
  File "realign_script.py", line 31, in <listcomp>
    [float(y) for y in x.split(' ')]
  File "realign_script.py", line 31, in <listcomp>
    [float(y) for y in x.split(' ')]
ValueError: could not convert string to float: ''

Using a slightly edited version of your script:

from cryosparc.tools import CryoSPARC
import json
import numpy as np
from pathlib import Path
from scipy.spatial.transform import Rotation as R
with open(Path('~/cs_instance_info.json').expanduser(), 'r') as f:
    instance_info = json.load(f)

cs = CryoSPARC(**instance_info)
assert cs.test_connection()

project_number = "P51"
workspace_number = "W6"

# src should be a symmetry expansion in C1 if you want
# to use the sym_expand/idx later
src_job_number = "J253"
src_job_particles_name = "particles"
src_job_vol_name = "volume"

lane_to_use = "default"

project = cs.find_project(project_number)
workspace = project.find_workspace(workspace_number)
job = project.find_job(src_job_number)
unrotated_particles = job.load_output(src_job_particles_name)

def parse_chimerax_matrix_string(chimx_mat:str) -> np.array:
    full_matrix = np.array(
        [
            [float(y) for y in x.split(' ')]
            for x in chimx_mat.split('\n')
        ]
    )
    rot_matrix = full_matrix[:,:3]
    translation = full_matrix[:, 3].T
    return (rot_matrix, translation)


chimx_mat = """   0.22976117   0.94489315   0.23320967 -41.36825376
    -0.95758188   0.17665770   0.22765985  29.48019271
     0.17391594  -0.27562475   0.94540163   5.72622007"""
rot_mat, translation = parse_chimerax_matrix_string(chimx_mat)

# we want a new center in voxels, not a translation in A
translation /= unrotated_particles['alignments3D/psize_A'][0]
new_box_center = unrotated_particles['blob/shape'][0][0] / 2 - translation

euler = R.from_matrix(rot_mat).as_euler('ZYZ')

# chimeraX translates *after* rotating, CryoSPARC translates *before*.
# thus, we need two VAT jobs. One rotating, one translating.
unshifted_job = workspace.create_job(
    type = "volume_alignment_tools",
    connections = {
        "particles": (src_job_number, src_job_particles_name),
        "volume": (src_job_number, src_job_vol_name)
    },
    params = {
        # yes, this is the right parameter name for the Euler angles
        'recenter_axang': ','.join(str(x) for x in euler)
    }
)
unshifted_job.queue(lane_to_use)
unshifted_job.wait_for_done()

unshift_id = unshifted_job.doc['uid']
shifted_job = workspace.create_job(
    type = "volume_alignment_tools",
    connections = {
        "particles": (unshift_id, "particles"),
        "volume": (unshift_id, "volume")
    },
    params = {
        "recenter_shift": ",".join(str(x) for x in new_box_center)
    }
)
shifted_job.queue(lane_to_use)
shifted_job.wait_for_done()

shifted_id = shifted_job.doc['uid']
shifted_results = project.find_job(shifted_id)
shifted_particles = shifted_results.load_output('particles')
shifted_particles['sym_expand/idx'] = 1
project.save_external_result(
    workspace_uid=workspace_number,
    dataset=shifted_particles,
    type='particle',
    slots=['sym_expand'],
    passthrough=(shifted_id, 'particles'),
    title='non-point-group sym expand'
)

Any idea what might be going wrong?

EDIT:

Adding an operation to strip trailing whitespace seemed to fix the issue, although the resulting volume alignment tools job does not have the expected rotation:

from cryosparc.tools import CryoSPARC
import json
import numpy as np
from pathlib import Path
from scipy.spatial.transform import Rotation as R
with open(Path('~/cs_instance_info.json').expanduser(), 'r') as f:
    instance_info = json.load(f)

cs = CryoSPARC(**instance_info)
assert cs.test_connection()

project_number = "P51"
workspace_number = "W6"

# src should be a symmetry expansion in C1 if you want
# to use the sym_expand/idx later
src_job_number = "J253"
src_job_particles_name = "particles"
src_job_vol_name = "volume"

lane_to_use = "default"

project = cs.find_project(project_number)
workspace = project.find_workspace(workspace_number)
job = project.find_job(src_job_number)
unrotated_particles = job.load_output(src_job_particles_name)

#def parse_chimerax_matrix_string(chimx_mat:str) -> np.array:
#    full_matrix = np.array(
#        [
#            [float(y) for y in x.split(' ')]
#            for x in chimx_mat.split('\n')
#        ]
#    )
#    rot_matrix = full_matrix[:,:3]
#    translation = full_matrix[:, 3].T
#    return (rot_matrix, translation)

def parse_chimerax_matrix_string(chimx_mat: str) -> np.array:
    print("Input string:")
    print(chimx_mat)

    lines = chimx_mat.strip().split('\n')  # Strip leading/trailing whitespace
    print("Split lines:")
    print(lines)

    full_matrix = np.array(
        [
            [float(y) for y in x.split()]  # No need to specify ' ' as split argument, it will split by any whitespace
            for x in lines
        ]
    )
    print("Parsed matrix:")
    print(full_matrix)

    rot_matrix = full_matrix[:, :3]
    translation = full_matrix[:, 3].T
    return rot_matrix, translation

chimx_mat = """ 0.22976117   0.94489315   0.23320967 -41.36825376
    -0.95758188   0.17665770   0.22765985  29.48019271
     0.17391594  -0.27562475   0.94540163   5.72622007 """
print(chimx_mat)
rot_mat, translation = parse_chimerax_matrix_string(chimx_mat)

# we want a new center in voxels, not a translation in A
translation /= unrotated_particles['alignments3D/psize_A'][0]
new_box_center = unrotated_particles['blob/shape'][0][0] / 2 - translation

euler = R.from_matrix(rot_mat).as_euler('ZYZ')

# chimeraX translates *after* rotating, CryoSPARC translates *before*.
# thus, we need two VAT jobs. One rotating, one translating.
unshifted_job = workspace.create_job(
    type = "volume_alignment_tools",
    connections = {
        "particles": (src_job_number, src_job_particles_name),
        "volume": (src_job_number, src_job_vol_name)
    },
    params = {
        # yes, this is the right parameter name for the Euler angles
        'recenter_axang': ','.join(str(x) for x in euler)
    }
)
unshifted_job.queue(lane_to_use)
unshifted_job.wait_for_done()

unshift_id = unshifted_job.doc['uid']
shifted_job = workspace.create_job(
    type = "volume_alignment_tools",
    connections = {
        "particles": (unshift_id, "particles"),
        "volume": (unshift_id, "volume")
    },
    params = {
        "recenter_shift": ",".join(str(x) for x in new_box_center)
    }
)
shifted_job.queue(lane_to_use)
shifted_job.wait_for_done()

shifted_id = shifted_job.doc['uid']
shifted_results = project.find_job(shifted_id)
shifted_particles = shifted_results.load_output('particles')
shifted_particles['sym_expand/idx'] = 1
project.save_external_result(
    workspace_uid=workspace_number,
    dataset=shifted_particles,
    type='particle',
    slots=['sym_expand'],
    passthrough=(shifted_id, 'particles'),
    title='non-point-group sym expand'
)

EDIT:

Ok got it to work! The input matrix has to be derived from measure rotation as applied to two duplicate maps - not models fit to the maps. Then, if the origin of both maps has been set to the center of the box, it does the trick! Thanks heaps!

1 Like
8

Meant to say - I was able to apply this approach to achieve local refinement of both monomers in 7PTY (EMPIAR-10828) and it gave a much nicer map (I suspect because it effectively reduces the orientational bias, as the monomers are asymmetrically disposed with respect to one another

2 Likes
9

Awesome! It’s always nice to hear that CryoSPARC tools enables these advanced workflows, and that they improve the final result! Glad we could help get it working for you :smile:

1 Like
10

Yes indeed! Hoping to use it in a workshop next month to show other folks how they can apply these kinds of workflows to their own samples :slight_smile:

4 Likes