Skip to content

Apply image shifting by fragment #124

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.

Sign up for GitHub

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

Open
hannahbaumann wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Apply image shifting by fragment #124

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
498c71f
Apply image shifting by fragment
hannahbaumann Jul 8, 2026
3fd980b
Update tests
hannahbaumann Jul 8, 2026
bba1661
Merge branch 'main' into align_fragment_based
hannahbaumann Jul 15, 2026
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
2 changes: 1 addition & 1 deletion src/openfe_analysis/rmsd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -328,12 +328,12 @@ def gather_rms_data(
skip = max(n_frames // 500, 1)

u_top = mda.Universe(pdb_topology)
u_top.select_atoms("protein").guess_bonds()

@IAlibay IAlibay Jul 16, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe you have to avoid doing this if there is no protein, otherwise it'll crash out.

Suggested change
u_top.select_atoms("protein").guess_bonds()
protein = u_top.select_atoms("protein")
if protein:
protein.guess_bonds()

Note: this also shows that we don't have any tests that are solvent only - we should raise an issue to deal with that.


for state_idx in range(n_lambda):
# cheeky, but we can read the PDB topology once and reuse per universe
# this then only hits the PDB file once for all replicas
universe = create_universe_single_state(u_top._topology, ds, state_idx)

prot = universe.select_atoms(protein_selection)
ligand = universe.select_atoms(ligand_selection)

Expand Down
3 changes: 3 additions & 0 deletions src/openfe_analysis/tests/test_rmsd_classes.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ def ligand(hybrid_system_skipped_pdb, simulation_skipped_nc):
universe = universe_utils.create_universe_single_state(
hybrid_system_skipped_pdb, simulation_skipped_nc, state=0
)
universe.select_atoms("protein").guess_bonds()
prot = universe.select_atoms("protein and name CA")
ligand = universe.select_atoms("resname UNK")
apply_transformations.apply_complex_alignment_transformations(universe, prot, [ligand])
Expand Down Expand Up @@ -150,6 +151,7 @@ def test_separate_ligands_fixes_pbc_spike(

# Combined approach should produce a spike
u_combined = universe_utils.create_universe_single_state(d["pdb"], ds, state=state_idx)
u_combined.select_atoms("protein").guess_bonds()
prot = u_combined.select_atoms("protein and name CA")
lig_A = u_combined.atoms[d["ligand_A_indices"]]
lig_B = u_combined.atoms[d["ligand_B_indices"]]
Expand All @@ -164,6 +166,7 @@ def test_separate_ligands_fixes_pbc_spike(

# Separate approach should fix it
u_separate = universe_utils.create_universe_single_state(d["pdb"], ds, state=state_idx)
u_separate.select_atoms("protein").guess_bonds()
prot = u_separate.select_atoms("protein and name CA")
lig_A = u_separate.atoms[d["ligand_A_indices"]]
lig_B = u_separate.atoms[d["ligand_B_indices"]]
Expand Down
3 changes: 3 additions & 0 deletions src/openfe_analysis/tests/utils/test_apply_transformations.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ def test_multichain_rmsd_shifting(simulation_skipped_nc, hybrid_system_skipped_p
u = universe_utils.create_universe_single_state(
hybrid_system_skipped_pdb, simulation_skipped_nc, 0
)
u.select_atoms("protein").guess_bonds()
prot = u.select_atoms("protein and name CA")
# Do other transformations, but no shifting
unwrap_tr = unwrap(prot)
Expand All @@ -75,7 +76,9 @@ def test_multichain_rmsd_shifting(simulation_skipped_nc, hybrid_system_skipped_p
u2 = universe_utils.create_universe_single_state(
hybrid_system_skipped_pdb, simulation_skipped_nc, 0
)
u2.select_atoms("protein").guess_bonds()
prot2 = u2.select_atoms("protein and name CA")
assert len(list(prot2.fragments)) == 2
apply_transformations.apply_complex_alignment_transformations(u2, protein=prot2)

R2 = rms.RMSD(prot2)
Expand Down
12 changes: 8 additions & 4 deletions src/openfe_analysis/utils/apply_transformations.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,11 +58,15 @@ def apply_complex_alignment_transformations(
# 1. Make molecules whole (protein + optional ligand)
transforms = [unwrap(group)]

@hannahbaumann hannahbaumann Jul 9, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This now makes the alignment very very slow. Unwrap loops over the fragments and before (when there was no bond information), there were a lot of fragments that were very small.

@IAlibay IAlibay Jul 16, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How much of a problem is this in practice?


# 2. Closest image shift for protein chains + ligand (if present)
chains = [seg.atoms for seg in protein.segments]
shift_targets = chains[1:] + ligands
# 2. Closest image shift for protein fragments + ligand (if present)
fragments = list(protein.fragments)

@IAlibay IAlibay Jul 16, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it worth checking for bonds in protein before you do this?

# Pick the largest fragment as the reference
ref_idx = max(range(len(fragments)), key=lambda i: fragments[i].n_atoms)
reference = fragments[ref_idx]
shift_targets = [f for i, f in enumerate(fragments) if i != ref_idx]
shift_targets += ligands
if shift_targets:
transforms.append(ClosestImageShift(chains[0], shift_targets))
transforms.append(ClosestImageShift(reference=reference, targets=shift_targets))

# 3. Align on protein backbone/atoms
transforms.append(Aligner(protein))
Expand Down
Loading