Parse hexadecimal resid from OpenMM in PDBParser

[rladbtls1998]

Is your feature request related to a problem?

OpenMM generate a PDB file with hexadecimal resid if there are more than 9999 residues.
However, in MDAnalysis, hexadecimal resids (treated as faulty resid) are set to 1, which causes problems in some cases.

Describe the solution you'd like

I'd like to change PDBParser._parseatoms to handle Hexadecimal resid, if there are alphabets in resid.
Changing

mdanalysis/package/MDAnalysis/topology/PDBParser.py

Lines 302 to 304 in c4af00b

	else:
	resid = int(line[22:26])
	# Wrapping

to

                    else:
                        if any([a.isalpha() for a in line[22:26]]):
                            resid = int(line[22:26], base=16) - 30960
                        else:
                            resid = int(line[22:26])
                        # Wrapping

Describe alternatives you've considered

Or setting unique resids (increment from last resid / preset failed resid starting number e.g. 10000) for residues with failed resid parsing can be another option (while this can be problematic for other cases).
Or editing Universe to get optional Parser(TopologyReaderBase) so one can give custom topology parser if needed.

Additional context

openmm.app.pdbfile.py : 484 (from OpenMM version 8.2.0)

def _formatIndex(index, places):
    """Create a string representation of an atom or residue index.  If the value is larger than can fit
    in the available space, switch to hex.
    """
    if index < 10**places:
        format = f'%{places}d'
        return format % index
    format = f'%{places}X'
    shiftedIndex = (index - 10**places + 10*16**(places-1)) % (16**places)  <----- This part
    return format % shiftedIndex

One major problem is that OpenMM uses decimal until resid 9999 and starts from hex A000 (== decimal 40960) in place of decimal resid 10000, so parser needs to subtract 30960 or just use resid 40960.
I think both may cause problems in other parts.

[aniketpati1121]

i can solve this issue

[rladbtls1998]

@aniketpati1121 while the issue itself is simple, I think it needs reviews or decisions from contributors.

[orbeckst]

If you look at #1897 and PR #1978 then I'd say MDAnalysis already supports hybrid36 PDB format.

Or is what you're looking for not "hybrid36"?

[rladbtls1998]

Thanks. didn't know about hybrid36 format.
However, OpenMM seems to use hexadecimal (base16), not hybrid36 format (base 36). (openmm/openmm#3501 and openmm/openmm#3522)
(in openmm/openmm#659 mentioned in #1897, they ended up implementing hex, not hybrid36)

I don't think it is practical to deal with every non-standard pdb format, so imo it is good enough to parse hex resid with hybrid36, since it prevents duplicate resids.
or
Letting Universe __init__ to accept optional parser (TopologyReaderBase) argument, so that users can provide custom parsers for non-standard formats, could also solve this issue.

Another minor problem is that hybrid36 PDB format is only used in atom serial number, but not in residue ids.
Simply changing

mdanalysis/package/MDAnalysis/topology/PDBParser.py

Lines 302 to 304 in c4af00b

	else:
	resid = int(line[22:26])
	# Wrapping

to

                else:
                   resid = hy36decode(4, line[22:26])
                    # Wrapping

does fix the missing hy36 resid issue.

[orbeckst]

Letting Universe __init__ to accept optional parser (TopologyReaderBase) argument, so that users can provide custom parsers for non-standard formats, could also solve this issue.

You can always provide a parser-subclass of TopologyReaderBase as a kwarg for Universe for topology_format=YourParser (and similar for the trajectory format=YourReader).

[orbeckst]

If I read your issue and the linked openmm issue/PR correctly then you're really asking for supporting reading of VMD-style hex numbered PDBs, right?

I don't think it is practical to deal with every non-standard pdb format, so imo it is good enough to parse hex resid with hybrid36, since it prevents duplicate resids.

Your suggestion would create another format, I think, because the resids would then have different numbers in MDAnalysis compared to how they were created.

If at all, I'd rather start from a defined standard. This might mean to analyze carefully what VMD and OpenMM are doing, write it up, and then decide if this means to support another parser. (We already have XPDB for extended residue numbers, we could have VMDPDB.)

[IAlibay]

Adding more support for new variants of an already deprecated file format just feels like spending efforts in the wrong direction.

Would time be better spent getting PDBx support into MDAnalysis? I know that OpenMM supports writing to PDBx.

[rladbtls1998]

Letting Universe __init__ to accept optional parser (TopologyReaderBase) argument, so that users can provide custom parsers for non-standard formats, could also solve this issue.

You can always provide a parser-subclass of TopologyReaderBase as a kwarg for Universe for topology_format=YourParser (and similar for the trajectory format=YourReader).

Thank you. I will do this.

If I read your issue and the linked openmm issue/PR correctly then you're really asking for supporting reading of VMD-style hex numbered PDBs, right?

I meant just using existing hy36 format implementation for resids.
Using currently existing hybrid36 to decode hex resids (from OpenMM/VMD etc) prevent duplicate resids by not failing and does not require any implementations except editing resid to also use hybrid36 .
(for context, in my case, all water hex resids failed to parse and were treated as resid 1, which caused problems due to duplicate resids)

Would time be better spent getting PDBx support into MDAnalysis? I know that OpenMM supports writing to PDBx.

There seems to be OpenMM converter that uses OpenMM to convert PDBx into MDAnalysis topology.

Adding more support for new variants of an already deprecated file format just feels like spending efforts in the wrong direction.

I agree.

FYI, after reading #1897 again, I noticed that #1897 was initially about hexadecimal format, not hybrid36. (similar to openmm/openmm#659 (comment))

Here's a really funny one packmol writes PDB indices in hexadecimal when they get too high. PDB really is the format that never stops giving

But I guess it never caused problems, since hy36 does decode hex into unique integers. (I'm not saying that hex format needs to be implemented though.)

TL;DR, I suggest changing

mdanalysis/package/MDAnalysis/topology/PDBParser.py

Lines 302 to 304 in c4af00b

	else:
	resid = int(line[22:26])
	# Wrapping

to use hybrid36 as I mentioned

[orbeckst]

Thanks for digging through all these old discussions. Amongst them was the comment #1897 (comment)

Does it make sense to allow this format also for residue ids?

The purist in me says "no, it's not a well-defined standard" ... but then I re-read the discussion on #1897 where we essentially did exactly what you're asking to do here: deal with odd files and just make them work.

I wholeheartedly agree that we should support PDBx/mmcif #2367 .

I also admit that your proposed fix seems quite straightforward and may not really influence much the speed at which issue #2367 is being solved. It also does not look as if it will break anything except that the actual resids will be numerically different (in base 10) from what VMD/openmm originally wrote.

In the light of the history of #1897 I am tending to support your change but I'd like to hear other @MDAnalysis/coredevs to weigh in — a simple 👍 or 👎 would be good, comments even better!

[BradyAJohnston]

My only concern would be in potential performance hit for parsing a large topology (likely minor), but ultimately the topology parsing is rarely the performance bottleneck so I think it's a sensible change.

[IAlibay]

I'm going to take a negative stance here. I understand that it's just a "small fix", but it's always "just another small fix" with PDB files, for some off standard adjustment that someone somewhere decided to randomly make because it was convenient for their tool at the time.

I will change my opinion to a positive one if:
0. Performance isn't affected for standard complying PDB parsing

1. This is the last off standard feature we allow in the PDB parser. If someone comes tomorrow with another feature that sits outside PDB V3, we reject it.
2. This isn't just a fix, but also a ton of documentation on what this new off standard feature is, where it comes from, and how it behaves.