Hi @rposert,
Hmm…I’m having trouble reproducing this. If I run a modern VAT job on particles from an old Homogeneous Refinement my outputs look as expected (like your first image). Could you make a new topic for this problem?
Split out as a separate thread here, thanks!
No, I never found a good dataset to test it on locally…but looking back at the code, could you try changing this line:
pose_recenter_vec_A = pose.apply(recenter_vector_A)
to
pose_recenter_vec_A = pose.inv().apply(recenter_vector_A)
This change unfortunately makes the resolution worse. I have a dataset I can provide for you to test with if useful? Where using the script to modify a defocus of a subparticle on the Z axis improves resolution substantially, but degrades resolution for a subparticle in the X-Y plane off-axis?
EDIT:
I added in some debugging and it seems like after combining a couple of the alterations in this thread, it is now giving max and min z-displacement values which match the measured distance in the structure . Running a test now on a case which should see improvement from this - will report back once the test is complete.
EDIT2:
That did the trick finally!! It is working now, even for off-axis subparticles! Improved resolution by 0.2Å for an off-axis subparticle! Here is the script I am using, thanks @rposert so much for your patience & perserverance!
from cryosparc.tools import CryoSPARC
import json
import numpy as np
from pathlib import Path
import re
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 = "P35"
workspace_number = "W1"
vol_align_uid = "J1081"
project = cs.find_project(project_number)
job = project.find_job(vol_align_uid)
particles = job.load_output("particles")
volume = job.load_output("volume")
apix = volume["map/psize_A"][0]
print("apix=",apix)
map_size_px = volume["map/shape"][0]
print("map_size_px=",map_size_px)
map_size_A = map_size_px * apix
if "recenter_shift" in job.doc["params_spec"]:
recenter_pattern = r"([\.\d]+), ?([\.\d]+), ?([\.\d]+) ?(px|A)?"
recenter_string = job.doc["params_spec"]["recenter_shift"]["value"]
print("recenter_string=",recenter_string)
recenter_regex = re.match(recenter_pattern, recenter_string)
print("len(recenter_regex.groups())=",len(recenter_regex.groups()))
print("group 4=",recenter_regex.group(4))
new_center_A = np.array([
float(x) for x in [
recenter_regex.group(1),
recenter_regex.group(2),
recenter_regex.group(3)
]
])
if recenter_regex.group(4) != "A":
new_center_A *= apix
print("new_center_if=",new_center_A)
print("new_center_A=",new_center_A)
else:
# rather than re-calculate the center of mass, just get it from the job log
# NOTE: cs.cli.get_job_streamlog() is hidden and experimental --
# it may change without warning
vol_align_log = cs.cli.get_job_streamlog(project_number, vol_align_uid)
centering_line = [x["text"] for x in vol_align_log if "New center will be" in x["text"]]
print("centering_line=",centering_line)
new_center_vox = re.search(r"([\.\d]+), ([\.\d]+), ([\.\d]+)", centering_line[0])
new_center_vox = np.array([
float(x) for x in [
new_center_vox.group(1),
new_center_vox.group(2),
new_center_vox.group(3)
]
])
new_center_A = new_center_vox * apix
# assume volume was well-centered
recenter_vector_A = new_center_A - map_size_A / 2
print("recenter_vector_A=",recenter_vector_A)
# recenter_vector_A is the vector pointing from the old center to the
# new center in the reference pose. Now we rotate this vector for each
# recentered particle by applying the pose
# iterating over rows is slow...
z_displacements = []
for row in particles.rows():
pose = R.from_rotvec(row["alignments3D/pose"])
pose_recenter_vec_A = pose.inv().apply(recenter_vector_A)
row["ctf/df1_A"] -= pose_recenter_vec_A[2]
row["ctf/df2_A"] -= pose_recenter_vec_A[2]
z_displacements.append(pose_recenter_vec_A[2])
max_z = np.max(z_displacements)
min_z = np.min(z_displacements)
mean_z = np.mean(z_displacements)
print(f"Max Z displacement: {max_z:.2f} Å")
print(f"Min Z displacement: {min_z:.2f} Å")
print(f"Mean Z displacement: {mean_z:.2f} Å")
out_job = project.create_external_job(
workspace_uid = workspace_number,
title = "Defocus-adjusted particles"
)
out_job.add_input(
type = "particle",
name = "input_particles"
)
out_job.connect(
target_input = "input_particles",
source_job_uid = vol_align_uid,
source_output = "particles"
)
out_job.add_output(
type = "particle",
name = "particles",
passthrough = "input_particles",
slots = ["ctf"],
alloc = particles
)
out_job.save_output("particles", particles)
out_job.log("Job created with adjust-defocus.ipynb")
out_job.stop()
I have for now left a few print statements etc in there for debugging purposes. @rposert it would be super helpful if this could be incorporated as an option in VAT!
EDIT3:
… although now the modified script no longer works for subparticles from the same assembly located on the Z axis, which did work with the original script, so there is still something wrong. Hmm. So now I have one script that works perfectly for subparticles in the X-Y plane, and one that works perfectly for subparticles that are on the Z-axis, but nothing that works for both. Clearly missing something…
EDIT4:
Ah I may have figured out what was going on. The subparticle in the X-Y plane was generated from a volume that had a different handedness from the one that generated the subparticle on the Z-axis. Testing to confirm that is the issue now and will report back.