Fix residue naming#2042
Conversation
|
pre-commit.ci autofix |
for more information, see https://pre-commit.ci
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #2042 +/- ##
==========================================
- Coverage 95.01% 90.70% -4.31%
==========================================
Files 212 213 +1
Lines 20572 20728 +156
==========================================
- Hits 19546 18801 -745
- Misses 1026 1927 +901
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
pre-commit.ci autofix |
for more information, see https://pre-commit.ci
IAlibay
left a comment
There was a problem hiding this comment.
This looks good thanks! - couple of things to fix but otherwise should be fine.
After this PR, will you be able to do the same for PlainMD, AFE, and SepTop?
|
|
||
| try: | ||
| data = rmsd.gather_rms_data(pdb_file, trj_file) | ||
| data = rmsd.gather_rms_data(pdb_file, trj_file, ligand_selection="resname LIG") |
There was a problem hiding this comment.
This will fail is someone provides their own residue names ahead of time (which is something we should encourage - i.e. power users should be rewarded not punished).
There are two potential solutions here:
- (lowest cost / temporary fix) We just store the alchemical residue names in the output of the setup unit and then fetch them in the analysis unit.
- (higher cost / more complicated) We pass in the alchemical components & the comp resids through all the way to the analysis unit and then extract the relevant residue names from that.
For now, I would suggest we just pass the residue names and have the residue name kwarg be optional and/or have a default of LIG, that way we don't break backwards compatibility.
We can clean up after the fact.
There was a problem hiding this comment.
Note that for hybrid toplogy, it's technically possible for both endstate molecules to have different names, so we should check that it works properly.
There was a problem hiding this comment.
I think a 3rd option might be best, where we do away with tracking the residue names and just keep the alchemical atom indices around. This would mean moving away from the gather_rms_data function and using the new openfe-analaysis API but we plan on doing that anyway, so this could be a good time to do it?
This would fix the case pointed out above and the fact that if a user sets the alchemical ligands and the cofactors to have the same residue name (LIG), we will still have the same issue in the analysis unit.
There was a problem hiding this comment.
I wouldn't consider the above blocking though we could have this as a temp fix and move to indices when we change the analysis for this protocol?
There was a problem hiding this comment.
Yes, I like the 3rd option as the long term fix, for now I went ahead with option 1, but happy to change it to a longer term fix.
Re two ligands different residue names: I tested this out and it assigns the residue name from the stateA ligand to the hybrid topology (in the HTF). We could still check for both residue names since it wouldn't hurt, but may also confuse the user?
There was a problem hiding this comment.
I think a 3rd option might be best, where we do away with tracking the residue names and just keep the alchemical atom indices around.
This somewhat fits within the longer term idea of switching comp_resids to just comp_indices (where you track the indices of each component in the system). However, that's a lot bigger lift.
I don't mind just tracking the alchemical indices in the setup results - that's something we nearly already do in the AFE Protocol. But yeah, for the sake of incremental PRs = faster velocity, I would vote for the quick fix now, then the medium/long term fix after.
Co-authored-by: Irfan Alibay <IAlibay@users.noreply.github.com>
|
Not blocking: Does a user have a way of setting the residue_name via a CLI input? It seems like the atom metadata fields are not preserved on writing an openff molecule to file? |
|
pre-commit.ci autofix |
for more information, see https://pre-commit.ci
I had forgotten that this was still an issue... The problem is that residue name is a per atom PDBResidueInfo property in rdkit molecules, and a) gufe SMCs don't serialize that, and b) it doesn't get written out when you write out an SDF file. See: OpenFreeEnergy/gufe#327 |
|
jthorton
left a comment
There was a problem hiding this comment.
This looks great @hannahbaumann though the role of LG2 is not clear currently as the output seems to use LG1 as the name for the hybrid ligand in the output topology. Having the other name around if we write out seperate PDBs for the end states would make sense but maybe its confusing until then?
|
|
||
| logger = logging.getLogger(__name__) | ||
|
|
||
| DEFAULT_LIGAND_RESNAMES = ("LG1", "LG2") |
There was a problem hiding this comment.
Sorry I'm going to go counter to what Josh said - I'm really not a big fan of this in such a big file. I'd prefer it being redefined than to have a hanging global variable that gets referenced effectively one time.
At worst let's just remove the default from run.
There was a problem hiding this comment.
Sorry @hannahbaumann I should have been more specific I did mean just declare in execute and remove the default from run if you agree to only have it defined in one place.
There was a problem hiding this comment.
Ok, removed the default from run now!
There was a problem hiding this comment.
Ah, but this now results in an API break...
There was a problem hiding this comment.
I think the API break is ok because we shouldn't have execution that spans multiple versions of openfe, but maybe I'm forgetting something!
|
pre-commit.ci autofix |
for more information, see https://pre-commit.ci
IAlibay
left a comment
There was a problem hiding this comment.
Sorry, I kinda convinced myself to go back to LIG 😅 - happy to discuss in a call if needed.
|
|
||
| if len(selection_indices) > 0: | ||
| sub_top = hybrid_topology.subset(selection_indices) | ||
| # Renumber sequentially so same-named residues (e.g. two COF cofactors, |
There was a problem hiding this comment.
Do we know if resids from the original PDB make it through all the way here?
I.e. if a user has a PDB file that starts at 4 instead of 1, until you do this renumbering - would the first residue also start at 4?
The main reason I ask this is that renumbering essentially loses any input information, so in some ways it makes the user experience worse because they can't then map things back to the input.
There was a problem hiding this comment.
Do the small molecules have a different chain ID from the protein residues? If so, would assigning the residue number at the same time we assign residue names to the offmol be suitable?
There was a problem hiding this comment.
I checked it and the small molecules are on a separate chain from the protein, so I'll switch to your suggestion and assign the residue number at the same time as the residue name!
| name = base if base not in used else _unique(used) | ||
| used.add(name) | ||
| else: | ||
| name = "COF" |
There was a problem hiding this comment.
[nit] this technically allows a case where the ligand(s) are named COF and the cofactors also get called COF.
Co-authored-by: Irfan Alibay <IAlibay@users.noreply.github.com>
|
pre-commit.ci autofix |
for more information, see https://pre-commit.ci
IAlibay
left a comment
There was a problem hiding this comment.
Looking great!
Couple of small things - mostly the test and the silent faillure with the get.
| name = _get_offmol_resname(small_mols[comp]) | ||
| assert name is not None | ||
| names.add(name) | ||
| alchem_resnames = sorted(names) |
There was a problem hiding this comment.
[nit] not something to fix now, but just an FYI that long term we might need to know what name belongs to what end state.
| pdb_file=pdb_file, | ||
| trajectory=trajectory, | ||
| checkpoint=checkpoint, | ||
| ligand_resnames=setup_results.outputs.get("alchemical_resnames"), |
There was a problem hiding this comment.
This will silently get None if the entry isn't there. it will eventually just make structural analysis fail completely silentlly (because it's alll wrapped in a try/except).
Co-authored-by: Irfan Alibay <IAlibay@users.noreply.github.com>
…y/openfe into fix_residue_naming
|
🚨 API breaking changes detected! 🚨 Griffe output |
Checklist
newsentry, or the changes are not user-facing.pre-commit.ci autofix.Manual Tests: these are slow so don't need to be run every commit, only before merging and when relevant changes are made (generally at reviewer-discretion).
Developers certificate of origin