Skip to content

[WIP] transformations#26

Open
JenkeScheen wants to merge 21 commits into
mainfrom
feat/transformations-first-steps
Open

[WIP] transformations#26
JenkeScheen wants to merge 21 commits into
mainfrom
feat/transformations-first-steps

Conversation

@JenkeScheen

@JenkeScheen JenkeScheen commented Jul 1, 2026

Copy link
Copy Markdown
Collaborator

Given selected Structures (one per ligand query), do the following transformations:

  • pick a reference protein conformation to be used from here on (best practices: largest ligand volume)
  • prepare protein (fix missing residues, heavy atoms, loops etc)
  • protonate protein and ligand (propka and dimorphite-DL, resp.)

side-track: also adding the full n=16 TYK2 set from the OFE industry benchmarks

@JenkeScheen
JenkeScheen force-pushed the feat/transformations-first-steps branch from 03f8080 to 4d23e35 Compare July 1, 2026 11:36
@JenkeScheen

Copy link
Copy Markdown
Collaborator Author

another side-track (sorry): while running TYK2 through the pipeline again, I ran into a pose that RDKit could not restore bond orders for because the connectivity was too wonky. I added a validator that catches this in a7324a3

@JenkeScheen

Copy link
Copy Markdown
Collaborator Author

@ianmkenney ready for a first look, in the meantime I'm just running TYK2 with this newest pipeline to see if it all works out.

@JenkeScheen
JenkeScheen requested a review from ianmkenney July 2, 2026 13:25
@JenkeScheen

Copy link
Copy Markdown
Collaborator Author

note that given a ligand series, there might be multiple net charges. Traditionally you would solve this by:

  • split the ligand series into a group per net charge
  • run an AlchemicalNetwork for each group
  • estimate DDGs for each group
  • run MLE for each group
  • merge group DGs back together

But this is a bit finnicky to set up robustly. OFE supports charge-changing transformations now, so we should just allow these for now -- it saves us from having to do the split/merge. We will need to 1) warn the user 2) make sure that in mapper there are enough edges between different charge state sections of the resulting network.

Comment thread src/openplaceholder/impl/transformations.py Outdated
Comment thread pyproject.toml Outdated
forcefield: str = "AMBER"


class MaxVolumeSiteTransformation(Transformation):

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

There is a lot going on in this transformation. I would intuitively expect that this only alters the structures based on ligand volumes but it also alters charge state. We are able to stack transformations.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

addressed by the split + explicit single-Structure helpers in ef4f9f5

out = Path(tmp) / "merged.pdb"
merged.atoms.write(str(out))
data = base64.b64encode(out.read_bytes()).decode()
return Structure(

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

enforces that the structure is always PDB after transformation even if the original structure was something else.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

addressed in ef4f9f5

def _setup(self) -> None:
pass

def _transform(self, structures: list[Structure]) -> list[Structure]:

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Tracing the execution path here is tricky and make futures maintainability difficult. The function names don't capture their actual behavior and they share too much responsibility.

Specifically, _prepare and _protonate silently rely on the correct ordering of the input sequence, and _prepare only operates on a single structure despite accepting a list.

We should rewrite this to be more explicit and decoupled. Something like this would be cleaner:

def _transform(self, structures) -> list[Structure]:
    chosen_complex = self._select_complex(structures)
    other_structures = [s for s in structures if s != chosen_complex]

    fixed_protein = self._fix_protein_heavy_atoms(chosen_complex)
    protonated_protein = self._protonate_protein_context(fixed_protein)
    ...

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Additionally, we should try to minimize when we write to disk. Aside from performance benefits, using the io module will narrow the classes of exceptions we might run into, like running out of space in /tmp/. I think the only place we fully need disk IO would be with pdb2pqr.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

addressed in ef4f9f5

)


class TestMaxVolumeSiteTransformation(TestCase):

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Good start but will need more

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

added some more in ef4f9f5

@JenkeScheen

Copy link
Copy Markdown
Collaborator Author

propka turns out not to be a suitable direction because it only predicts pkas and there is no turnkey solution to actually using those predictions to protonate the whole complex in one go. Instead we'll go with what's been built in protonate-utils and leave the separate protein and ligand protonation methods as interchangeable modules. This does present a new problem though, because these are separate protonators there may be 'clashes' in the form of e.g. hbond donor-donor interactions between the ligand and the protein. We'll use prolif and some internal heuristics to find and fix them.

@JenkeScheen

Copy link
Copy Markdown
Collaborator Author

@ianmkenney can you have another look? this PR has become too inflated at this point imo

@ianmkenney ianmkenney left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I think we need to look again at how these transformations fit into the overall pipeline and how they should be stacked.

Given a pre-validated set of structures, we first apply MaxVolumeSiteSubstitutionTransformation (or whatever we would call it), where the "canonical structure" is immediately substituted into all other structures. From there we can apply the ProteinPreparationTransformation. By caching, we can avoid redoing the extra N-1 calculations. We then apply the ComplexProtonationTransformation. This completely removes the need to keep track of an implicit structure ordering dependency.

structures
)

self.assertEqual(len(structures), 1)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I'm not sure how much this is adding over individual transformation tests. You could remove the MaxVolumeSiteSelectionTransformation.transform call everything would still pass. Both of your final asserts are targeted to a single transformation but the way it's written doesn't attribute the change to either of the function calls.


selected = transformation.transform([small, large])

self.assertEqual(selected[0].ligand_name, "large")

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

why do we expect a reordering of structures after a transformation?

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

oh I see now...

atom.name,
)
continue
clashing_h = min(p_hs, key=lambda i: float(np.linalg.norm(p_xyz[i] - l_xyz[l_heavy])))

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

It seems to me that this could strip a neutral lysine down to a single hydrogen by mistake.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

These should probably be in the ComplexProtonationTransformation class as methods instead of individual classes.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I'd say hold off on doing this type of abstraction until we know we need it. As some of the other comments go over, most of this behavior should probably be packaged into the ComplexProtonationTransformation class. There may come a day when we add more implementations and something like this will be helpful.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

As stated in other comments, not sure this subpackage is needed.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I feel like this should all be part of the ComplexProtonationTransformation class.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Let's go over some of the logic here in our next sync

logger = logging.getLogger(__name__)


def _perceive_ligand(structure: Structure) -> Chem.Mol | None:

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

This is called in multiple validators. In a new PR I'll have this called in the parent class and cache the results.

opt into only the stages they want -- e.g. prepare without protonating, or
protonate a set that was selected elsewhere.

Ordering contract: **index 0 is the canonical protein context.** The selection

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

We should not use order dependence here.

out_file = pdb_io.PDBFile()
out_file.set_structure(structure)
out_pdb = Path(tmp) / "protein_h.pdb"
out_file.write(str(out_pdb))

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Can this be done in-memory?

Comment thread src/openplaceholder/impl/protonation/reconcile.py Outdated

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Don't do ordering, stick to duplication of effort. We can handle it in a later PR with caching.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants