Extending PARSE and AMBER FF for non-standard amino-acids

[stefdoerr]

Hi, I am working on extending the forcefields in pdb2pqr for non-standard amino-acids. I read my structure in with PARSE and output AMBER naming. I am not using pdb2pqr for creating PQR files, I just want to protonate and optimize my protein and standardize the atom names.

For example this protein contains a TYS amino-acid.

For now I've managed to roughly create the XML and DAT files for it (excluding charge/radius) TYS.zip which I will merge into AA.xml and PARSE.dat, however I have some questions

1. Do I need charges/radii for pdb2pqr to protonate my structures? Or are those necessary just for generating PQR files (which I don't need)?
2. If I do need them, where can I obtain the radii / charges from? Could I just take them from GAFF2 / AM1-BCC, would that work for PARSE?
3. I want my output FF to be AMBER, however since I am creating these residues myself I could use the same naming scheme for PARSE+AMBER. In that case do I still need to modify AMBER.DAT/names or would pdb2pqr keep the residue unmodified if it cannot find it in the AMBER ff?
4. Do I need anything in the .name files? As far as I can tell if I use the same naming scheme in PARSE and AMBER it should not be necessary?

Thanks in advance!

[stefdoerr]

On my first attempt I got the following error:

~/Work/pdb2pqr/pdb2pqr/biomolecule.py in create_residue(self, residue, resname)
    970                     klass = getattr(aa, resname)
    971                 except AttributeError:
--> 972                     klass = getattr(na, resname)
    973                 residue = klass(residue, refobj)
    974         except (KeyError, NameError):

AttributeError: module 'pdb2pqr.na' has no attribute 'TYS'

Does that mean I'd need to add a new class to aa.py for each AA I am adding to the forcefield? That makes it a bit hard to extend, unless I find some way of dynamically adding classes to the pdb2pqr.aa module?

Just to test, I added it like this and it worked

class TYS(Amino):
    """TYS class."""

    def __init__(self, atoms, ref):
        Amino.__init__(self, atoms, ref)
        self.reference = ref

    def letter_code(self):
        return "X"

Seems to protonate the residue fine, here you can see TYS together with the previous residue 362

[sobolevnrm]

You are correct; you current have to add a new TYS class to make this work. We're trying to refactor the code to use CIF containers as the base object which should make this easier. However, that work is ongoing and we don't have an ETA right now. @intendo

[sobolevnrm]

Here are some answers to your questions above:

1. Do I need charges/radii for pdb2pqr to protonate my structures? Or are those necessary just for generating PQR files (which I don't need)?

No, those aren't needed for protonation.

2. If I do need them, where can I obtain the radii / charges from? Could I just take them from GAFF2 / AM1-BCC, would that work for PARSE?

That seems like a good starting place.

3. I want my output FF to be AMBER, however since I am creating these residues myself I could use the same naming scheme for PARSE+AMBER. In that case do I still need to modify AMBER.DAT/names or would pdb2pqr keep the residue unmodified if it cannot find it in the AMBER ff?

I think that should work but am not 100% sure. That part of the code is very messy and it's been a long time since I've looked at it.

4. Do I need anything in the `.name` files? As far as I can tell if I use the same naming scheme in PARSE and AMBER it should not be necessary?

I agree: I think it is not necessary.

[stefdoerr]

Great! Thanks a lot for the help :) I'll try to make it work within the current limitations.

[shaoqx]

Thank you @stefdoerr for the amazing research on this! I just ran into this problem today with MSE.
Would it be great if we can submit PRs that add stuffs into these 3 files ( aa.py AA.xml PARSE.dat ) so that PDB2PQR will support more and more modified amino acid? (e.g.: enforcing using standard names from PDB)
I also agree it will be great if user can just support several files in the input.

[shaoqx]

For any future reader of this discussion, note that you need to align the coordinate with the standard backbone coordinates in pdb2pqr when making the XML file. Otherwise, the H on the backbone will fly away in the protonated structure.
This is a PDB file I made as the reference of the alignment. One needs to change the coordinate and rotate the $\psi$ $\phi$ angle.

HETATM    1  N   PSD P   1       1.201   0.847   0.000  0.00  0.00      PSDO N  
HETATM    2  CA  PSD P   1       0.000   0.000   0.000  0.00  0.00      PSDO C  
HETATM    3  C   PSD P   1      -1.250   0.881   0.000  0.00  0.00      PSDO C  
HETATM    4  O   PSD P   1      -2.185   0.660  -0.784  0.00  0.00      PSDO O  
HETATM    5  C-1 PSD P   1       2.339   0.216   0.000  0.00  0.00      PSDO C  
HETATM    6  N+1 PSD P   1      -1.252   1.877   1.023  0.00  0.00      PSDO N  
HETATM    7  H   PSD P   1       1.201   1.847   0.000  0.00  0.00      PSDO H