Secondary Structure Detection Fails Occasionally

[BradyAJohnston]

Secondary Structure Detection Fails Occasionally

[OLaprevote]

I already saw this with other structures. I think it's a limitation from biotite itself, more precisely the algorithm they use for secondary structure detection (P-SEA) that is using only the dihedral angles to detect SS, contrary to DSSP which uses hydrogen bonds as well.

I see two ways out of this:

• Go the hardcore open-source warrior way and implement DSSP into biotite. Technically there is already a python implementation of DSSP, actually containing two: calculations can be done with either numpy or pytorch. So the work would be mainly to port this to biotite. I would actually be willing to get this done as I think it would benefit the structural biology community as a whole in the long term, but I'd say if there is no hidden difficulty it's still a day work (+ doing the PR and getting it accepted).
• I think it was already suggested in an issue I can't find, but reading the secondary structure attribute out of the mmCIF file itself may be another way out. It's probably less tricky. Back then I chose not to go down this road so that any kind of file (+ eventually MD trajectories) could get the SSE, but, well, it can't hurt.

[BradyAJohnston]

Just had a go at implementing a purely Geometry Nodes implementation of DSSP and it actually works pretty well:

Currently it detects helices (alpha and 3-10 helices) by using their hydrogen bonding, then beta-strands based on their angles. This doesn't check at all for hydrogen bonding of the beta-sheets. I think that will be a limitation until 3.6 and more nodes come along. Seems to work 'well enough' on the structures that I have tested though. Given this is purely for visualisation, I'm less determined to get it 100% accurate.

I think I would still much prefer to get access to the author-deposited SS from the file, but this can potential work for MD simulations and when there isn't SS available.

[tubiana]

I found a bug but it's coming from BIOTITE example and not for b-sheets...
They have classified G and I as coil instead of helices.
see here : #203 (comment)

Cheers!
Thibault.

[OLaprevote]

Yup, just ran 3AEH_A into the sse tutorial of biotite and the problem definitely comes from biotite, maybe (probably?) from the SS detection algorithm 🙃
(I may be wrong and it's something happening in biotite itself, but I also saw that alpha helices were sometimes shorter than when using DSSP, so that's my main guess)

[BradyAJohnston]

We could also implement the SS detection algorithm inside of GeometryNodes itself. This might be slightly trickier, but if it can be done then each frame of an MD animation can have SS detection applied to it, or even if you deform a protein then you can update the SS.

[OLaprevote]

Implementing either P-SEA-like or DSSP-like algorithm in Blender Geometry Nodes would be… wow.

• P-SEA "only" requires to compute phi and psi. Given I could do a N->O vector it should be doable. I imagine defining the right SS phi/psi regions could be the difficult part, though phi/psi "squares" at first may work.
• DSSP would be a bigger chunk of math to implement, with some fun matricial operations (pyDSSP detection is ~100 lines and quite clean). Let's just say I'm not sure it's doable in one full day of work. Having a first look at the algorithm it looks feasible, or at the very least worth a try x)

It's not just dynamics: it would easily be computed for lots of different cases, e.g. in the long run maybe allow people to "draw" proteins with secondary structures.

For MD, though, I have to make a case: it should also be possible to implement pyDSSP in MDAnalysis. Same argument than for biotite: it would be useful not just for MolecularNodes but for a broader community, and you having a cool & precise user case is a good nudge for this to happen. I would gladly do it but it may take a few months (you saw my handling of the SS attribute here).
Getting DSSP working in MDAnalysis was already a requested feature but the only attempt to do so was getting a conda implementation of C++ DSSP code to work in there and it was deemed too impractical as it introduced way too many dependencies.

Anyways I'd say both projects are broader than just a quick bug fix. Do you think right now it would be easy to somehow allow to "paint" the protein backbone with SS attribute, so that users could quickly change a few amino-acids to the right SS when it's not to their convenience? (Which… could actually be really cool for lots of different applications)

[BradyAJohnston]

If possible, having DSSP (or something at least mostly approximating it) inside of geometry nodes would be great. Might not be possible until the Index of Nearest node comes in potentially 3.6.

Certainly in the short term, accessing the SS that is already stored in a file is likely the best option, rather than trying to re-compute it. I haven't looked into it for biotite or MDAnalysis to see if either readers give access to that information, but in theory it should be.

[BradyAJohnston]

I'm probably going to consider this solved for now as I implemented taking SS from the mmtf file when downloading from the PDB. The calculations will still fail a bit when opening a local file. The MOL_utils_dssp node will help in a lot of those cases and I will work on improving it also.