Hello, gromacs (and chemistry) newbie here. I have a set of NMR pdb structures where I am labelling each residue with its secondary structure. I initially used gmx DSSP for this as it’s the most widely used algorithm, however when plotting phi and psi angles coloured by assignment (using only a-helix and b-strand), it is just completely random. MDAnalysis’s DSSP does not exhibit this behaviour and can distinguish almost perfectly between the structures on the ramachandran plot, as well as delta 2D (which uses chemical shifts instead of pdb). I am confident that the code is correct between algorithms as they all have very similar output.
As taken from the gmx dssp documentation:
Note that gmx dssp currently is not capable of reproducing the secondary structure of proteins whose structure is determined by methods other than X-ray crystallography (structures in .pdb format with incorrect values in the CRYST1 line) due to the incorrect cell size in such structures.
I’ve looked at a few papers that the current version of gmx dssp is based on, but I think I am limited by my lack of knowledge of X-ray crystallography as my studies so far have been purely NMR-based. I’m aware that the algorithm was trained on crystal structures as per the DSSP 4 publication, but I’m confused as to why the CRYST1 values appear to have such a huge impact to render the output useless.
My question here is: how does the gmx dssp algorithm use the CRYST1 cell size in it’s SS determination, and why does the absence of correct CRYST1 values appear to completely mess up the calculation?
here is an example of the obvious mismatch (I have around 6,000 pdbs and they all have this profound disagreement)
I’ve just run my pdbs on DSSP v.4 (which the gmx version is based on), and I’m getting good accuracy. It is only the gromacs implementation that the CRYST1 values seems to be implemented in.
Having summed up all of the above, it is necessary to clarify that to get the same results as from the DSSP v.4 algorithm, it is mandatory to use the following parameters: “-hmode dssp -hbond energy -clear -pihelix”. Also, to get the same results as from the DSSP v.2 algorithm, the following parameters are needed: “-hmode dssp -hbond energy -clear -pihelix -nopolypro”.
But I have not found this to be the case. I am thoroughly confused, my pdbs are correctly formatted (with the default CRYST1 values), and I have played around with the different flag options for gmx and they all seem to produce completely random assignments.
I have not heard of any issues with gmx dssp. Is it not simply that there are incorrect values in the CRYST1 record? What happens when you remove this line?
Well… I’m now very weirded out because when I removed the CRYST1 line in a sample pdb, the agreement appeared to increase compared to DSSP.v4 - HOWEVER I have somehow lost 20 residues…
gmx dssp (CRYST1 line still there): (length 105)
~BTTTSS~BS~SSB~~S~~S~~S~~S~~SSS~S~SS~S~S~~~~SSSB~SS~S~TTTTSSSS~S~~~~S~B~S~SS~~SS~SSSSS~~~S~TTT~SBB~BTTB~~
gmx dssp CRYST1 line removed: (length 85)
~PSTTEEEEEESSS~EEEEEEPPSSSEEEEEEE~TTS~EEEEEE~SS~SEEEE~STT~EEEEEEEEEBTTBEEEEEE~SS~SS~~
DSSP.v4 (CRYST1 line still there): (length 105)
-PS--TT----EEEEEES--SS-EEEEEEPPSS--SEEEEEEE-TTS-EEEEEE-SS-SEEEE-S--TT-EEEEEEEEEBTTB----EEEEEE---S--S-SS--
I want to note that I’ve since run all three (MDAnalysis dssp, DSSP.v4, and gmx DSSP) on a subset of 11 crystal structures (ie, structures with accurate CRYST1 values) and the agreement was great. So this appears to be an issue in how gmx dssp handles the CRYST1 line that was changed from DSSP.v4 that has not been documented anywhere that I can find. But I am also hesitant to claim this as I’ve never seen anyone else have this issue with NMR structures.
I do not believe that the mismatch is caused by errors in the CRYST1 line as I have around 6,000 pdbs, all realised via NMR, so they all have a CRYST1 line that looks like: CRYST1 1.000 1.000 1.000 90.00 90.00 90.00 P 1 1
That probably does cause it. GROMACS will read this as a box of size 0.1 nm. We should probably interpret such CRYST1 records as no periodicity. Removing these lines should fix your issue.
Apologies for my late reply. I’ve tested removing CRYST1 and changing the box size to 1000 on three random pdbs. I’m honestly baffled - both of these changes caused the agreement to dramatically improve on the assignment of SS. I understand the logic that removing the CRYST1 line might improve agreement as the absence of box size may cause the algorithm to activate whatever default box is optimised to perform best (a complete guess, obviously not 0.1x0.1x0.1nm), but I do not understand why an absurdly large box size also increases agreement (unless there is a max size where the default box is applied if exceeded - again a total guess).
I’ve rerun gmx dssp on the untouched pdbs so many times and yes, the disagreement shown below perseveres. I do think that my other comment where a bunch of residues disappeared was probably an error on my part somewhere as I haven’t been able to recreate it, so apologies.
original agreement for pdb 1DXW.pdb:
d2d: -----S----SEEEEE-SS-EEEEEEETTEEEEEHHHHSS-----SSS--PPSEEEGGGTEEE---STT---BB----S--EEEE--SS---EEEEEEETTEEEESS-B-GGGGTTTTT--EE-TTSBEEEEEEEEEESSSSEEEEEEEEHHHHHHHHH-S----
mda: -----------EEEEE----EEEEEEE--EEEEEHHHH---------------EEEHHH-EEE---------EE-------EEEE-------EEEEE-E--EEEE---E-HHHH-------EE----EE-EEE-EEEE----E-EEEEEEHHHHHHHHH------??
dssp.v4: -----S----SEEEEE-SS-EEEEEEETTEEEEEHHHHSS-----SSS--PPSEEEGGGTEEE---STT---BB----S--EEEE--SS---EEEEEEETTEEEESS-B-GGGGTTTTT--EE-TTSBEEEEEEEEEESSSSEEEEEEEEHHHHHHHHH-S----
gmx: ~BTTB~SSSS~SSSS~S~SS~SSSS~S~SS~~STTTSB~TTSSS~SSSBB~~~BSB~SBTTS~SS~SSSSBSSSSSSSSS~SSS~SSSSBSSSS~BSSS~~S~S~S~S~S~STTTTSSS~SSSTT~BS~SSSSSSBSSSBSS~SSSBBSSS~S~SSSSSSS~~~~
gmx dssp with CRYST1 line removed:
~~~~~S~~~~SEEEEE~SS~EEEEEEETTEEEEEHHHHSS~~~~~SSS~~PPSEEEGGGTEEE~~~STT~~~BB~~~~S~~EEEE~~SS~~~EEEEEEETTEEEESS~B~GGGGTTTTT~~EE~TTSBEEEEEEEEEESSSSEEEEEEEEHHHHHHHHH~S~~~~
gmx dssp with CRYST1 line set to 1000:
~~~~~S~~~~SEEEEE~SS~EEEEEEETTEEEEEHHHHSS~~~~~SSS~~PPSEEEGGGTEEE~~~STT~~~BB~~~~S~~EEEE~~SS~~~EEEEEEETTEEEESS~B~GGGGTTTTT~~EE~TTSBEEEEEEEEEESSSSEEEEEEEEHHHHHHHHH~S~~~~
original agreement for pdb 2NBV.pdb:
d2d: --SEEEEEEEEEEETTEEEEP-SPEEEEEEE-TT-PEEEEEEETTT--EEEEEEE-TTSEEEEEPSS-SSS-EEEEEESSSS-EEEEEE--S-GGGHHHHHHHHHHHHH--------EEEEEEE-SS-EEEEEEETT--HHHHHHHHHHHT-S-SSS-EEE-SBSEESTTS-IIIII-SSSEEEEEE--
mda: ---EEEEEEEEEEE--EEEE----EEEEEEE-----EEEEEEE-----E--EEEE----EEEEE--------EEEEEE-----EEEEEE----HHHHHHHHHHHHHHHH-------EEEEEEE----EEEEEEE----HHHHHHHHHHH--------EEE--E-EE-----HHHHH----EEEEEE--
dssp.v4: --SEEEEEEEEEEETTEEEEP-SPEEEEEEE-TT-PEEEEEEETTT--EEEEEEE-TTSEEEEEPSS-SSS-EEEEEESSSS-EEEEEE--S-GGGHHHHHHHHHHHHH--------EEEEEEE-SS-EEEEEEETT--HHHHHHHHHHHT-S-SSS-EEE-SBSEESTTS-IIIII-SSSEEEEEE--
gmx: ~BBTSSS~BSS~SSS~SSBTTTBSSS~S~SSSS~SS~~SS~~~SS~SS~SSS~S~~SSSSS~S~SSBS~SSSBSSS~~S~TTTTTSSSSSSSSS~S~SSSSSSS~S~S~~S~SSSSSS~SSSS~STTTSSSS~SSSSSSSSSS~SSSS~SSSSSS~~S~~SSSS~S~SS~SSSS~S~SS~SBS~SS~~
gmx dssp with CRYST1 line removed:
~~SEEEEEEEEEEETTEEEEP~SPEEEEEEE~TT~PEEEEEEETTT~~EEEEEEE~TTSEEEEEPSS~SSS~EEEEEESSSS~EEEEEE~~S~GGGHHHHHHHHHHHHH~=~~~~~~EEEEEEE~SS~EEEEEEETT~~HHHHHHHHHHHT~S~SSS~EEE~SBSEESTTS~IIIII~SSSEEEEEE~~
gmx dssp with CRYST1 line set to 1000:
~~SEEEEEEEEEEETTEEEEP~SPEEEEEEE~TT~PEEEEEEETTT~~EEEEEEE~TTSEEEEEPSS~SSS~EEEEEESSSS~EEEEEE~~S~GGGHHHHHHHHHHHHH~=~~~~~~EEEEEEE~SS~EEEEEEETT~~HHHHHHHHHHHT~S~SSS~EEE~SBSEESTTS~IIIII~SSSEEEEEE~~
original agreement for pdb 6Z0G.pdb:
d2d: --GGGT--SS-HHHHHHHHHHHHHHHHHH-HHHHHHHH-SS-----
mda: --HHH------HHHHHHHHHHHHHHHHHH-HHHHHHHH--------
dssp.v4: --GGGT--SS-HHHHHHHHHHHHHHHHHH-HHHHHHHH-SS-----
gmx: ~HHHHTSSSSS~SSSSS~SSSTT~S~SSSS~S~~~S~~BB~S~B~~
gmx dssp with CRYST1 line removed:
~~GGGT~~SS~HHHHHHHHHHHHHHHHHH~HHHHHHHH~SS~~~~~
gmx dssp with CRYST1 line set to 1000:
~~GGGT~~SS~HHHHHHHHHHHHHHHHHH~HHHHHHHH~SS~~~~~
If it truly is the case that GROMACS does read and take the CRYST1 line into account in its prediction (which would have been an undocumented change from DSSP.v4), I think it makes sense to request an amendment where if the box size is exactly the default of 1, 1, 1, gmx dssp ignores it (or treats it as below the minimum size, much like 1000,1000,1000 exceeds the max size). This would allow for the algorithm to be used for NMR based analysis.
I don’t know of any change in behavior. All GROMACS tools apply periodic boundary conditions (by default), as otherwise many results would be wrong when analyzing simulation trajectories.
Using a CRYST1 of 1,1,1 to indicate no periodic boundary conditions is a very nasty hack (not surprising of course, knowing the PDB format and its history). As it seems very unlikely that people would use an 0,1 nm large box, we can flag such cases a no PBC when reading PDB.
PS: GROMACS uses a zero size box for indicating no PBC when we do not have a PDB flag. That is a much better solution than using 1 …