InSty - changes in calcStatisticsInteractions

prody/proteins/interactions.py

@@ -46,7 +46,7 @@
            'calcHydrogenBondsTrajectory', 'calcHydrophobicOverlapingAreas',
            'Interactions', 'InteractionsTrajectory', 'LigandInteractionsTrajectory',
            'calcSminaBindingAffinity', 'calcSminaPerAtomInteractions', 'calcSminaTermValues',
-           'showSminaTermValues', 'showPairEnergy']
+           'showSminaTermValues', 'showPairEnergy', 'checkNonstandardResidues']
 
 
 def cleanNumbers(listContacts):
@@ -160,6 +160,53 @@ def get_energy(pair, source):
     import numpy as np
     import importlib.resources as pkg_resources    
 
+    aa_correction = {
+        # Histidine (His)
+        'HSD': 'HIS',   # Protonated at ND1 (HID in AMBER)
+        'HSE': 'HIS',   # Protonated at NE2 (HIE in AMBER)
+        'HSP': 'HIS',   # Doubly protonated (HIP in AMBER)
+        'HID': 'HIS',   # AMBER name, protonated at ND1
+        'HIE': 'HIS',   # AMBER name, protonated at NE2
+        'HIP': 'HIS',   # AMBER name, doubly protonated
+        'HISD': 'HIS',  # GROMACS: protonated at ND1
+        'HISE': 'HIS',  # GROMACS: protonated at NE2
+        'HISP': 'HIS',  # GROMACS: doubly protonated
+
+        # Cysteine (Cys)
+        'CYX': 'CYS',   # Cystine (disulfide bridge)
+        'CYM': 'CYS',   # Deprotonated cysteine, anion (GROMACS)
+
+        # Aspartic acid (Asp)
+        'ASH': 'ASP',   # Deprotonated Asp (GROMACS: ASH, AMBER: AS4)
+
+        # Glutamic acid (Glu)
+        'GLH': 'GLU',   # Deprotonated Glu (GROMACS: GLH, AMBER: GL4)
+
+        # Lysine (Lys)
+        'LYN': 'LYS',   # Deprotonated lysine (GROMACS: LYN, AMBER: LYP)
+
+        # Arginine (Arg)
+        'ARN': 'ARG',   # Deprotonated arginine (rare, GROMACS)
+
+        # Tyrosine (Tyr)
+        'TYM': 'TYR',   # Deprotonated tyrosine (GROMACS: TYM)
+
+        # Serine (Ser)
+        'SEP': 'SER',   # Phosphorylated serine (GROMACS/AMBER)
+
+        # Threonine (Thr)
+        'TPO': 'THR',   # Phosphorylated threonine (GROMACS/AMBER)
+
+        # Tyrosine (Tyr)
+        'PTR': 'TYR',   # Phosphorylated tyrosine (GROMACS/AMBER)
+
+        # Non-standard names for aspartic and glutamic acids in low pH environments
+        'ASH': 'ASP',   # Protonated Asp
+        'GLH': 'GLU',   # Protonated Glu
+    }
+
+    pair = [aa_correction.get(aa, aa) for aa in pair]    
+
     try:
         # Python 3
         with pkg_resources.path('prody.proteins', 'tabulated_energies.txt') as file_path:
@@ -187,11 +234,45 @@ def get_energy(pair, source):
     return data_results
 
 
+def checkNonstandardResidues(atoms):
+    """Check whether the atomic structure contain non-standard residues and inform to replace the name 
+    to the standard one to be that non-standard residues are treated in a correct way while computing 
+    interactions.
+
+    :arg atoms: an Atomic object from which residues are selected
+    :type atoms: :class:`.Atomic`
+    """
+
+    try:
+        coords = (atoms._getCoords() if hasattr(atoms, '_getCoords') else
+                    atoms.getCoords())
+    except AttributeError:
+        try:
+            checkCoords(coords)
+        except TypeError:
+            raise TypeError('coords must be an object '
+                            'with `getCoords` method')
+
+    amino_acids = ["ALA", "ARG", "ASN", "ASP", "CYS", "GLU", "GLN", "GLY", "HIS", "ILE", 
+                   "LEU", "LYS", "MET", "PHE", "PRO", "SER", "THR", "TRP", "TYR", "VAL"]
+
+    aa_list = atoms.select('name CA').getResnames()
+    nonstandard = []
+
+    for i in aa_list:
+        if i not in amino_acids:
+            nonstandard.append(i)
+
+    LOGGER.info('There are several non-standard residues in the structure.')
+    LOGGER.info('Replace the non-standard name in the PDB file with the equivalent name from the standard one if you want to include them in the interactions.')
+    LOGGER.info("Residues: {0}.".format(' '.join(nonstandard)))   
+
+
 def showPairEnergy(data, **kwargs):
     """Return energies when a list of interactions is given. Energies will be added to each pair of residues 
     at the last position in the list. Energy is based on the residue types and not on the distances.
     The unit of energy is kcal/mol. The energies defined as 'IB_nosolv', 'IB_solv' are taken from [OK98]_ and 
-    'CS' from InSty paper.
+    'CS' from InSty paper (under preparation). 
 
     :arg data: list with interactions from calcHydrogenBonds() or other types
     :type data: list
@@ -2007,6 +2088,7 @@ def calcStatisticsInteractions(data, **kwargs):
                     "energy": get_energy([element.split('-')[0][:3], element.split('-')[1][:3]], energy_list_type)
                 }
             except:
+                LOGGER.warn('energy information is not available for ', element.split('-')[0][:3], element.split('-')[1][:3])
                 stats[element] = {
                     "stddev": np.round(np.std(values),6),
                     "mean": np.round(np.mean(values),6),

prody/proteins/tabulated_energies.txt

@@ -398,120 +398,3 @@ PRO LYS 0.3	-0.5	0.0
 PRO ARG -0.17	-2.43	0.0
 PRO HIS 0.33	-1.8	0.0
 PRO PRO -0.03	-0.83	0.0
-GLY HSE 0.01	-2.47	0.0
-ALA HSE 0.01	-3.29	0.0
-VAL HSE -0.23	-4.14	-1.1
-ILE HSE -0.02	-4.55	-1.1
-LEU HSE 0.25	-4.85	-1.1
-CYS HSE -0.64	-4.79	-0.8
-MET HSE -0.17	-4.47	-0.5
-PHE HSE 0.79	-3.82	-0.6
-TYR HSE 0.34	-3.78	-1.1
-TRP HSE -0.05	-4.5	-1.7
-SER HSE -0.38	-3.12	-0.5
-THR HSE -0.05	-2.73	-0.5
-ASP HSE 0.0	-2.93	-0.4
-ASN HSE -0.52	-3.05	-1.2
-GLU HSE -0.1	-3.15	-0.4
-GLN HSE -0.31	-4.2	-1.2
-LYS HSE -0.01	-2.14	0.0
-ARG HSE 0.35	-3.24	-0.4
-HSE GLY 0.01	-2.47	0.0
-HSE ALA 0.01	-3.29	0.0
-HSE VAL -0.23	-4.14	-1.1
-HSE ILE -0.02	-4.55	-1.1
-HSE LEU 0.25	-4.85	-1.1
-HSE CYS -0.64	-4.79	-0.8
-HSE MET -0.17	-4.47	-0.5
-HSE PHE 0.79	-3.82	-0.6
-HSE TYR 0.34	-3.78	-1.1
-HSE TRP -0.05	-4.5	-1.7
-HSE SER -0.38	-3.12	-0.5
-HSE THR -0.05	-2.73	-0.5
-HSE ASP 0.0	-2.93	-0.4
-HSE ASN -0.52	-3.05	-1.2
-HSE GLU -0.1	-3.15	-0.4
-HSE GLN -0.31	-4.2	-1.2
-HSE LYS -0.01	-2.14	0.0
-HSE ARG 0.35	-3.24	-0.4
-HSE HSE 0.38	-3.08	-0.5
-HSE PRO 0.33	-1.8	0.0
-PRO HSE 0.33	-1.8	0.0
-GLY HSD 0.01	-2.47	0.0
-ALA HSD 0.01	-3.29	0.0
-VAL HSD -0.23	-4.14	-1.1
-ILE HSD -0.02	-4.55	-1.1
-LEU HSD 0.25	-4.85	-1.1
-CYS HSD -0.64	-4.79	-0.8
-MET HSD -0.17	-4.47	-0.5
-PHE HSD 0.79	-3.82	-0.6
-TYR HSD 0.34	-3.78	-1.1
-TRP HSD -0.05	-4.5	-1.7
-SER HSD -0.38	-3.12	-0.5
-THR HSD -0.05	-2.73	-0.5
-ASP HSD 0.0	-2.93	-0.4
-ASN HSD -0.52	-3.05	-1.2
-GLU HSD -0.1	-3.15	-0.4
-GLN HSD -0.31	-4.2	-1.2
-LYS HSD -0.01	-2.14	0.0
-ARG HSD 0.35	-3.24	-0.4
-HSD GLY 0.01	-2.47	0.0
-HSD ALA 0.01	-3.29	0.0
-HSD VAL -0.23	-4.14	-1.1
-HSD ILE -0.02	-4.55	-1.1
-HSD LEU 0.25	-4.85	-1.1
-HSD CYS -0.64	-4.79	-0.8
-HSD MET -0.17	-4.47	-0.5
-HSD PHE 0.79	-3.82	-0.6
-HSD TYR 0.34	-3.78	-1.1
-HSD TRP -0.05	-4.5	-1.7
-HSD SER -0.38	-3.12	-0.5
-HSD THR -0.05	-2.73	-0.5
-HSD ASP 0.0	-2.93	-0.4
-HSD ASN -0.52	-3.05	-1.2
-HSD GLU -0.1	-3.15	-0.4
-HSD GLN -0.31	-4.2	-1.2
-HSD LYS -0.01	-2.14	0.0
-HSD ARG 0.35	-3.24	-0.4
-HSD HSD 0.38	-3.08	-0.5
-HSD PRO 0.33	-1.8	0.0
-PRO HSD 0.33	-1.8	0.0
-GLY HSP 0.01	-2.47	0.0
-ALA HSP 0.01	-3.29	0.0
-VAL HSP -0.23	-4.14	-1.1
-ILE HSP -0.02	-4.55	-1.1
-LEU HSP 0.25	-4.85	-1.1
-CYS HSP -0.64	-4.79	-0.8
-MET HSP -0.17	-4.47	-0.5
-PHE HSP 0.79	-3.82	-0.6
-TYR HSP 0.34	-3.78	-1.1
-TRP HSP -0.05	-4.5	-1.7
-SER HSP -0.38	-3.12	-0.5
-THR HSP -0.05	-2.73	-0.5
-ASP HSP 0.0	-2.93	-0.4
-ASN HSP -0.52	-3.05	-1.2
-GLU HSP -0.1	-3.15	-0.4
-GLN HSP -0.31	-4.2	-1.2
-LYS HSP -0.01	-2.14	0.0
-ARG HSP 0.35	-3.24	-0.4
-HSP GLY 0.01	-2.47	0.0
-HSP ALA 0.01	-3.29	0.0
-HSP VAL -0.23	-4.14	-1.1
-HSP ILE -0.02	-4.55	-1.1
-HSP LEU 0.25	-4.85	-1.1
-HSP CYS -0.64	-4.79	-0.8
-HSP MET -0.17	-4.47	-0.5
-HSP PHE 0.79	-3.82	-0.6
-HSP TYR 0.34	-3.78	-1.1
-HSP TRP -0.05	-4.5	-1.7
-HSP SER -0.38	-3.12	-0.5
-HSP THR -0.05	-2.73	-0.5
-HSP ASP 0.0	-2.93	-0.4
-HSP ASN -0.52	-3.05	-1.2
-HSP GLU -0.1	-3.15	-0.4
-HSP GLN -0.31	-4.2	-1.2
-HSP LYS -0.01	-2.14	0.0
-HSP ARG 0.35	-3.24	-0.4
-HSP HSP 0.38	-3.08	-0.5
-HSP PRO 0.33	-1.8	0.0
-PRO HSP 0.33	-1.8	0.0