First version of coordinate generatign processor #328
Conversation
There was a problem hiding this comment.
Overall looks pretty good in terms of a draft. In my opinion the main points for making it a working part of vermouth are:
- missing docstrings and cryptic variable names. It is very time consuming to review (especially) the math if it is not clearly defined what even the shape of a variable is.
- more clear signature of functions; in several cases functions return and do multiple things depending on the Type of the input variables and/or combination of input variables. I think in several places this can be more clear.
- more testing. All 400 lines of code will need both unit tests and integral tests. This will be a huge amount of work, so probably this should be starting soon.
- bench marking and profiling; there are several functions that operate on the complete graph. It doesn't look to me like they are very expensive to evaluate, but to be honest martinize2 cannot really afford to become much slower. So I'd highly recommend profiling in terms that that processor scales linearly with graph length/number of edges at the very lest.
| def mindist(A, B): | ||
| return ((A[:, None] - B[None])**2).sum(axis=2).min(axis=1) |
There was a problem hiding this comment.
I understand that for many linear algebra functions it is tempting to use "math like" names. However, it is very hard to read. I'd suggest to at least in the doc string and even better the variable name itself, to at least mention if the objects A,B are both vectors, matrices and what the requirements on their shapes are. Ideally you'd even use 'vector_a' or 'matrix_a'. This holds throughout the document.
There was a problem hiding this comment.
Personally I really disagree that math like names and algorithms are harder to read, for sure in (almost) purely mathematical algorithms. I have more difficulty with convoluted prose. I do agree that the concise variables should be explicated in the docstring, along with the equations used. The same for mathematical texts, right? In general, especially for straightforward mathematical algorithms, I do prefer staying close to the equations, as these will also be present in accompanying documentation (including the article and docs) and overall it makes code more digestable. This does require a consistent notation, which is a bit harder in code, and it may be better to use simple prefixes (a for array, m for matrix, v for vector, etc). I'll update the docstrings. If you guys prefer prefixes, then I can add them. But longer names would not generally help, because lines become longer, may need to be split and in some cases may cause more use of memory and/or overhead.
There was a problem hiding this comment.
I'm somewhere in the middle. Changing A, B to vector_A, vector_B indeed does not add much over describing them in the docstring. However, A and B (usually) refer to some physical thing, such as an angle or position, and adding those does benefit understanding/comprehension. So how about e.g. position_A or pos_A, etc?
But longer names would not generally help, because lines become longer, may need to be split and in some cases may cause more use of memory and/or overhead.
I don't mind long lines too much (anything under 100 characters is fine, over 120 is not) (comment and docstrings should be <80). And I don't think your memory/overhead argument holds. In this function the intermediate results of the subtraction, sum, and min are all stored in memory separately at some point. Maybe also the squaring even.
There was a problem hiding this comment.
In the function mentioned they are fully abstract; they're just sets of vectors, 2D numpy arrays, matrices. They could be positions, angles, RNA transcript counts, whatever. In the functions _chiral and _out, the input signatures refer to anchor and target, and the arithmetic short-hand variables refer to the schematics. I think that warrants the use as is. The one-letter references to conic point sets and some other things could benefit from less cryptic variable names.
| Determine particles with known coordinates connected to | ||
| particles with given indices. |
| return P, S | ||
|
|
||
|
|
||
| def triple_cone(n=24, distance=0.125, angle=109.4): |
There was a problem hiding this comment.
Can't we sum up the cone functions such that you say cones(n_conscutive_SP3s, anchor ...)? And do I get this right that this is for linear bonding i.e. triple cone is SP3-SP3-SP3 not SP3-(SP3)2? So branching doesn't get covered?
There was a problem hiding this comment.
It will for sure be possible to generalize that. It's not, however, just for linear bonding. You can select multiple points (branches) at each level, which will give you (access to) the corresponding possible points at the next.
|
|
||
| b1 | ||
| \ | ||
| b2--a-->X | ||
| / | ||
| b3 |
There was a problem hiding this comment.
Would be nice if the variable names would at least be matching the drawing.
| if control == None: | ||
| B = [ | ||
| molecule.nodes[ix]['position'] | ||
| for ix, v in molecule[anchor].items() | ||
| if v != {} | ||
| ] | ||
| else: | ||
| B = [ molecule.nodes[ix]['position'] for ix in control ] |
There was a problem hiding this comment.
I suppose that the shape of B determines how the atom is rebuild?
Now the default is to get B from the anchor, which otherwise you get from indices. For making the code more readable and testable, I think it would be better to either accept anchors or only accept indices. It seems redundant to have both. Maybe I'm mistaken, but I'm happy to learn.
There was a problem hiding this comment.
Agreed. It's a remnant, but I was afraid of removing it in case it proved useful later.
| if target is not None: | ||
| molecule.nodes[target]['position'] = pos |
There was a problem hiding this comment.
Also here I'd argue it is better to either update the molecular position or to return the built position. That would be a much clearer signature.
There was a problem hiding this comment.
I think this too is a remnant of some earlier snippets. What would be preferred? I think changing by reference is maybe too implicit, and the function would be more useful if it only returns positions, but that adds more overhead in the calling function.
| if control is None: | ||
| B = [ | ||
| molecule.nodes[ix]['position'] | ||
| for ix, v in molecule[anchor].items() | ||
| if v != {} | ||
| ] | ||
| else: | ||
| B = [ molecule.nodes[ix]['position'] for ix in control ] |
There was a problem hiding this comment.
same as above; also the fact that this code fragment occurs in two functions tells me that it could also be recycled.
There was a problem hiding this comment.
@Tsjerk how thoroughly did you test this code on pdb files?
I tested it with several protein PDB files (single chain and multichain), including with listed missing atoms and with all side chains simply removed by hand. In all cases the atoms are rebuild to the point that conversion to CG is feasible. Some particles, notably hydrogens, may be misplaced, but that doesn't significantly affect the resulting CG model. It does, at present, make it impossible to use this to fix structures for atomistic simulations. I did not test cases with non-protein missing atoms.
There was a problem hiding this comment.
In fact the placement of hydrogen matters quite a lot in Martini3. For example, rings mapped from COG without hydrogen will be too small. In general I agree that is fine because energy minimization takes care of it, but the elastic network is a problem. The EN will refer to the "unrelated" coordinates. How sever is this problem? I don't know. I do know, however, that for Martin3 DNA it has gone to a point that martinize2 is unusable with the desired mapping, because the EN counteracts the relaxed volume geometry. Therefore accurate reconstruction of atoms does matter to some extent also for CG models.
Furthermore, martinize2 was written with the intention of being usable independent of resolution. Yet this coordinate generator is another thing demonstrating that it is actually targeted specifically for the conversion of AA to CG nothing more nothing less. And in this aspect it is even worse than PDB2GMX. Not that we would have time to fix all this, but it should have been clear that an approximate coordinate generator is not really good enough for what martinze2 claims to be. That is my opinion on it.
There was a problem hiding this comment.
Agreed, the method needs to be improved. This is a first take on it. It is, however, not an 'approximate coordinate generator'. To get this first version, a few shortcuts were taken, which can be relatively easily expanded on, but that takes time. The elastic network could be a concern. Is that applied on side chains? Or only on backbone beads? The latter would not be a concern. The first would, because the orientation of side chains cannot be considered optimized with respect to the neighborhood. A workaround would be to do a simulation with backbones fixed to allow side chains to equilibrate before applying the EN. But that's workflow-wise for the case where you have (many) missing parts). I wouldn't want to suggest doing too much energy minimizing stuff in martinize.
A somewhat greater concern of mine is the lack of chirality information, which is really problematic for atomistics scale models. But that cannot only be solved in this part, but also needs directives in the blocks/links that can then be taken into account. Here, the processor already has the mechanics (up/down), but nothing to base the placement on.
If it is only hydrogens that need to be built, the result is pretty much the same as what pdb2gmx does. pdb2gmx cannot build sidechains, not even a single atom. It can build you a C- or N-terminus, but that's it. So this is already miles better.
Okay, an inventory of the problems:
- The CB is, for some reason, considered SP2, which causes HB1 (HB in VAL) to be placed on top of CA. This only happens if CB is built by the routine, not if it is already present. (expected to be easy to fix)
- Rings are not constructed correctly. (expected to be somewhat hard)
- The arginine guanidinium group is not planar. (expected to be relatively easy)
Otherwise the sidechains are well-formed, although not always placed optimally. There are optimizations possible without going to energy minimization.
There was a problem hiding this comment.
This sounds pretty good to me.
I hereby veto any and all energy minimisation in vermouth/martinize2, which should answer the fist paragraph. It would, however, be good (excellent even) to issue a warning when reconstructing chains of atoms (a H is fine, a CH3 probably also fine, CH2-CH3 probably not) to advise users to energy minimize before doing e.g. elastic networks.
Chirality: Detect when something could be chiral (95%), and warn. No longer our problem. We can revisit this later.
1 and 3 should be fixed.
For 2 we can cheat as far as I'm concerned. Just assume it's a perfect polygon, and reconstruct the ring as-a-whole. Not ideal, especially regarding sugar, but maybe good enough?
Bottom line: best effort, and if we/you don't trust the result make a loud warning explaining to the user what they can do. If there's a warning no files are written (unless maxwarn), so it'll no longer be our problem.
| try: | ||
| rup, rdown = distance | ||
| except TypeError: | ||
| rup = rdown = distance |
There was a problem hiding this comment.
in terms of signature this isn't very clear either, because once distance acts as a float as in the default but then it can also be a list/array? You could have it as tuple or dict or any other way as well. In addition this try accept seems very broad and prevents sensible errors if the Type of distances is actually wrong.
There was a problem hiding this comment.
A distance could be specified now as float or as two floats. There's a laziness there. So for signatures sake it shall be two floats always.
| if up is not None: | ||
| molecule.nodes[up]['position'] = pup | ||
| if down is not None: | ||
| molecule.nodes[down]['position'] = pdown | ||
|
|
||
| return pup, pdown |
There was a problem hiding this comment.
Same as before; either return distances or update molecule.
| '''Class for subgraphs of (missing) atoms with anchors''' | ||
| def __init__(self, molecule, limbs): | ||
| self.molecule = molecule | ||
| self.anchors = limbs[0][0] | ||
| self.bases = limbs[0][1] |
There was a problem hiding this comment.
Especially, needs doc-string of what is what
| def get_missing_atoms(molecule): | ||
| ''' | ||
| Determine particles without coordinates in molecule. | ||
| ''' | ||
| return [ | ||
| ix for ix in molecule | ||
| if molecule.nodes[ix].get('position') is None | ||
| ] |
There was a problem hiding this comment.
Replace for use of selector_has_position https://github.com/marrink-lab/vermouth-martinize/blob/master/vermouth/selectors.py#L67
| ''' | ||
| return { | ||
| a for ix in indices for a in molecule[ix].keys() | ||
| if molecule.nodes[a].get('position') is not None |
There was a problem hiding this comment.
Add select_missing to selector.py?
There was a problem hiding this comment.
Is an option, but isn't easier to just negate the outcome of has_position?
There was a problem hiding this comment.
here yes, but in the other function the purpose is to select all atoms without coordinates. The purpose suggests that selectors.py is a better place. Probably less so for 'anchors' which are more specific to the building routines.
There was a problem hiding this comment.
Makes sense, especially in combination with selectors.filter_minimal. In general I'm rather unhappy with the selection functions at the moment, since it's completely unclear whether they select molecules, atoms, nodes in molecules, or something else entirely. So feel free to use your own best judgement, and we'll clean up the selectors at a later moment.
| if coords is not None: | ||
| P[mindist(P, coords) < contact**2] = np.nan | ||
| P = P.reshape((len(self.limbs), -1, 3)) | ||
| for px in P.transpose((1,0,2)): |
|
|
||
| prev = len(missing) | ||
| missing = get_missing_atoms(molecule) | ||
| #print(len(missing)) |
| return molecule | ||
|
|
||
|
|
||
| def main(args): |
There was a problem hiding this comment.
Should be (re)moved before merging
There was a problem hiding this comment.
In general I like it!
I do agree with Fabian's comments though, and on top of that I have a few nitpicks. Also, it may be nice if this code produces some user output (logging!).
All that said, with deadlines and such approaching I'm open to merging this as-is (with a TODO comment, maybe) to speed up development at the cost of accruing maintenance debt
@Tsjerk how thoroughly did you test this code on pdb files?
@fgrunewald @jbarnoud Opinion?
As it is completely disconnected from the rest, I am not against it if it accelerates dev. Though, it means that actual integration will require more work since at least the tests are not in a usable state. |
|
Another thing that needs doing before merging is to add this to bin/martinize2 so it's actually used.
How so? There are no tests for this code, which means Travis is still happy |
|
I am open for the option of merging without the code-coverage being 100%. However, I think it should come with at least basic testing on the geometry functions and someone testing of the model in the workflow of martinize. As far as I understand it is not integrated in to bin/martini2 yet and has only been tested via the main module provided. Last time I tried running the main program it got tangled up in dependency errors. Maybe that is just me calling it; should be like this or am I mistaken: Running the basic processor to reconstruct H atoms on CH2 groups it also gave a bunch of strange coordinate choices. Probably when you map them it is fine but it looks super strange in of you visualize it. Bottom line can someone else try running some tests. It could be that my PDBs are little special but as of now I'd really recommend an additional layer of testing. |
That's the thing: if it get integrated to bin/martinize2, then it needs tests. Once it is in the command line, then it can mess up the user's science. Or... it comes behind a flag with a massive warning about it being under development. |
Sounds reasonable. I'm a little surprised it struggles with hydrogens, I would expect those to be the easiest case. Could you try reconstructing the hydrogens from a small protein, and seeing whether the resulting structure survives Gromacs minimization? If so, I'd say is perfectly fine enough for now. |
If it's only hydrogens, it works fine. That is indeed the simples case. There's a flaw with like the methyl of isoleucine.
These are built, but distorted, making them unsuitable for atomistic simulations, but doable for coarse graining. |
No thorough benchmark, but it takes 260 ms on my machine (md41) to rebuild 18273 atoms in a protein that is stripped to the bare backbone. |
|
@Tsjerk Any progress/update? |
There was a problem hiding this comment.
Much better already :)
I still don't like the single letter variables. For example in _chiral and _out you can rename a -> anchor (yes, override the argument) and B -> bases. But I won't die on this hill.
The docstrings you wrote I like (nit aside). It's quite important that at least all public methods (e.g. mindist) describe the arguments like you do with _chiral. For _private methods we can be a little bit more relaxed. But ideally all have detailed docstrings.
| } | ||
|
|
||
|
|
||
| def align_z(v): | ||
| ''' | ||
| Generate rotation matrix aligning z-axis onto v (and vice-versa). | ||
| Generate rotation matrix aligning z-axis onto vector v (and vice-versa). |
There was a problem hiding this comment.
Is it an idea to rename v to vector here? It's a bit of a nitpick from my side, admittedly
| # normalized vectors, the alignment is a 180 degree | ||
| # rotation around the resultant vector. |
There was a problem hiding this comment.
This comment is so good I wonder if it should be part of the docstring.
| Returns | ||
| ------- | ||
| None | ||
| position |
There was a problem hiding this comment.
Nit, also for the other docstring below:
Returns
----------
np.array
The position of the outward pointing atom.
Sphinx is rather picky
The coordinate generating processor. Unfortunately, the chirality is not always consistent, even for backbone, which needs to be fixed, and the side chains with rings are not well constructed. Otherwise, the coordinates are such that a protein stripped from side chains is rebuilt to be processable with martinize2.