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,17 +160,63 @@ 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:
-            data = np.loadtxt(file_path, skiprows=1, dtype=str)
+            data = np.loadtxt(file_path, dtype=str)
     except: 
         # Python 2.7
         import pkg_resources
         file_path = pkg_resources.resource_filename('prody.proteins', 'tabulated_energies.txt')
         with open(file_path) as f:
-            data = np.loadtxt(f, skiprows=1, dtype=str)
-
+            data = np.loadtxt(f, dtype=str)
 
     sources = ["IB_nosolv", "IB_solv", "CS"]
     aa_pairs = []
@@ -180,21 +226,66 @@ def get_energy(pair, source):
 
     lookup = pair[0]+pair[1]
 
-    return data[np.where(np.array(aa_pairs)==lookup)[0]][0][2:][np.where(np.array(sources)==source)][0]
+    try:
+        data_results = data[np.where(np.array(aa_pairs)==lookup)[0]][0][2:][np.where(np.array(sources)==source)][0]
+    except ImportError:
+        raise ImportError('Please replace non-standard names of residues with standard names.')
+
+    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 XX and 
-    'CS' from YY.
+    The unit of energy is kcal/mol. The energies defined as 'IB_nosolv', 'IB_solv' are taken from [OK98]_ and 
+    'CS' from InSty paper (under preparation). 
 
     :arg data: list with interactions from calcHydrogenBonds() or other types
     :type data: list
 
     :arg energy_list_type: name of the list with energies 
                             default is 'IB_solv'
     :type energy_list_type: 'IB_nosolv', 'IB_solv', 'CS'
+
+
+    .. [OK98] Keskin O., Bahar I., Badretdinov A.Y., Ptitsyn O.B., Jernigan R.L., 
+    Empirical solvet-mediated potentials hold for both intra-molecular and 
+    inter-molecular inter-residues interactions,
+    *Protein Science* **1998** 7: 2578–2586.
     """
 
     if not isinstance(data, list):
@@ -1988,12 +2079,21 @@ def calcStatisticsInteractions(data, **kwargs):
     for element in elements:
         if element not in stats:
             values = [t[1] for t in interactions_list if t[0] == element]
-            stats[element] = {
-                "stddev": np.round(np.std(values),6),
-                "mean": np.round(np.mean(values),6),
-                "weight": np.round(float(len(values))/len(data), 6),
-                "energy": get_energy([element.split('-')[0][:3], element.split('-')[1][:3]], energy_list_type)
-            }
+
+            try:
+                stats[element] = {
+                    "stddev": np.round(np.std(values),6),
+                    "mean": np.round(np.mean(values),6),
+                    "weight": np.round(float(len(values))/len(data), 6),
+                    "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),
+                    "weight": np.round(float(len(values))/len(data), 6)
+                }
 
     statistic = []
     for key, value in stats.items():
@@ -2002,8 +2102,11 @@ def calcStatisticsInteractions(data, **kwargs):
             LOGGER.info("  Average [Ang.]: {}".format(value['mean']))
             LOGGER.info("  Standard deviation [Ang.]: {0}".format(value['stddev']))
             LOGGER.info("  Weight: {0}".format(value['weight']))
-            LOGGER.info("  Energy [kcal/mol]: {0}".format(value['energy']))
-            statistic.append([key, value['weight'], value['mean'], value['stddev'], value['energy']])
+            try:
+                LOGGER.info("  Energy [kcal/mol]: {0}".format(value['energy']))
+                statistic.append([key, value['weight'], value['mean'], value['stddev'], value['energy']])
+            except:
+                statistic.append([key, value['weight'], value['mean'], value['stddev']])
         else: pass
 
     statistic.sort(key=lambda x: x[1], reverse=True)