Skip to content

Commit

Permalink
Merge pull request prody#2000 from jamesmkrieger/insty_energy
Browse files Browse the repository at this point in the history 
more Insty energy things
  • Loading branch information
@jamesmkrieger
jamesmkrieger authored Dec 1, 2024
2 parents 7d4ab54 + 3a1f6b1 commit 3e3d822
Show file tree
Hide file tree
Showing 3 changed files with 128 additions and 28 deletions.
36 changes: 35 additions & 1 deletion prody/dynamics/plotting.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@
'showPairDeformationDist','showMeanMechStiff',
'showPerturbResponse', 'showTree', 'showTree_networkx',
'showAtomicMatrix', 'pimshow', 'showAtomicLines', 'pplot',
'showDomainBar']
'showDomainBar', 'showSelectionMatrix']


def showEllipsoid(modes, onto=None, n_std=2, scale=1., *args, **kwargs):
Expand Down Expand Up @@ -2403,3 +2403,37 @@ def showTree_networkx(tree, node_size=20, node_color='red', node_shape='o',
showFigure()

return mpl.gca()

def showSelectionMatrix(matrix, atoms, selstr_x=None, selstr_y=None, **kwargs):
"""
Show a matrix similarly to showAtomicMatrix but only for
selected atoms based on *selstr_x* and *selstr_y*
:arg selstr_x: a selection string used with sliceAtomicData with axis=0
to slice the matrix in the x axis and label it with corresponding atoms
:type selstr_x: str
:arg selstr_y: a selection string used with sliceAtomicData with axis=1
to slice the matrix in the y axis and label it with corresponding atoms
:type selstr_y: str
If either of these are left as *None*, then no slicing is performed in that direction.
"""
atoms_x = atoms
atoms_y = atoms

if selstr_x is not None:
if not isinstance(selstr_x, str):
raise TypeError('selstr_x should be a str')

matrix = sliceAtomicData(matrix, atoms, selstr_x, axis=0)
_, atoms_x = sliceAtoms(atoms, selstr_x)

if selstr_y is not None:
if not isinstance(selstr_y, str):
raise TypeError('selstr_y should be a str')

matrix = sliceAtomicData(matrix, atoms, selstr_y, axis=1)
_, atoms_y = sliceAtoms(atoms, selstr_y)

return showAtomicMatrix(matrix, atoms=[atoms_x, atoms_y], **kwargs)
112 changes: 85 additions & 27 deletions prody/proteins/interactions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
from numpy import *
from prody import LOGGER, SETTINGS, PY3K
from prody.atomic import AtomGroup, Atom, Atomic, Selection, Select
from prody.atomic import flags
from prody.atomic import flags, sliceAtomicData
from prody.utilities import importLA, checkCoords, showFigure, getCoords
from prody.measure import calcDistance, calcAngle, calcCenter
from prody.measure.contacts import findNeighbors
Expand Down Expand Up @@ -277,9 +277,11 @@ def checkNonstandardResidues(atoms):
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' (non-solvent-mediated), 'IB_solv' (solvent-mediated)
are taken from [OK98]_ and 'CS' from InSty paper (under preparation).
The unit of energy is kcal/mol. The energies defined as 'IB_nosolv' (non-solvent-mediated) and
'IB_solv' (solvent-mediated) are taken from [OK98]_ and 'CS' from InSty paper (under preparation).
Protonation of residues is not distinguished. The protonation of residues is not distinguished.
'IB_solv' and 'IB_nosolv' have RT units and 'CS' has units of kcal/mol.
Known residues such as HSD, HSE, HIE, and HID (used in MD simulations) are treated as HIS.
:arg data: list with interactions from calcHydrogenBonds() or other types
Expand Down Expand Up @@ -2133,9 +2135,9 @@ def calcStatisticsInteractions(data, **kwargs):
(2) average distance of interactions for each pair [in Ang],
(3) standard deviation [Ang.],
(4) Energy [in kcal/mol] that is not distance dependent. Energy by default is solvent-mediated
from [OK98]_ ('IB_solv'). To use non-solvent-mediated entries ('IB_nosolv') from [OK98]_ or
solvent-mediated values obtained for InSty paper ('CS', under preparation) change
`energy_list_type` parameter.
from [OK98]_ ('IB_solv') in RT units. To use non-solvent-mediated (residue-mediated) entries ('IB_nosolv')
from [OK98]_ in RT units or solvent-mediated values obtained from MD for InSty paper ('CS', under preparation)
in kcal/mol, change `energy_list_type` parameter.
If energy information is not available, please check whether the pair of residues is listed in
the "tabulated_energies.txt" file, which is localized in the ProDy directory.
Expand Down Expand Up @@ -2192,14 +2194,15 @@ def calcStatisticsInteractions(data, **kwargs):
}

statistic = []
unit = 'RT' if energy_list_type in ['IB_solv', 'IB_nosolv'] else 'kcal/mol'
for key, value in stats.items():
if float(value['weight']) > weight_cutoff:
LOGGER.info("Statistics for {0}:".format(key))
LOGGER.info(" Average [Ang.]: {}".format(value['mean']))
LOGGER.info(" Standard deviation [Ang.]: {0}".format(value['stddev']))
LOGGER.info(" Weight: {0}".format(value['weight']))
try:
LOGGER.info(" Energy [kcal/mol]: {0}".format(value['energy']))
LOGGER.info(" Energy [{0}]: {1}".format(unit, value['energy']))
statistic.append([key, value['weight'], value['mean'], value['stddev'], value['energy']])
except:
statistic.append([key, value['weight'], value['mean'], value['stddev']])
Expand Down Expand Up @@ -3464,15 +3467,10 @@ def buildInteractionMatrixEnergy(self, **kwargs):
:type energy_list_type: str
'IB_solv' and 'IB_nosolv' are derived from empirical potentials from
O Keskin, I Bahar and colleagues from [OK98]_.
O Keskin, I Bahar and colleagues from [OK98]_ and have RT units.
'CS' is from MD simulations of amino acid pairs from Carlos Simmerling
and Gary Wu.
.. [OK98] Keskin O, Bahar I, Badretdinov AY, Ptitsyn OB, Jernigan RL,
Empirical solvent-mediated potentials hold for both intra-molecular
and inter-molecular inter-residue interactions
*Protein Sci* **1998** 7(12):2578-2586.
and Gary Wu in the InSty paper (under preparation) and have units of kcal/mol.
"""

import numpy as np
Expand Down Expand Up @@ -3791,7 +3789,10 @@ def showCumulativeInteractionTypes(self, **kwargs):
:arg DiBs: score per disulfide bond
:type DiBs: int, float
:arg selstr: selection string for focusing the plot
:type selection: str
:arg energy: sum of the energy between residues
default is False
:type energy: bool
Expand All @@ -3805,16 +3806,25 @@ def showCumulativeInteractionTypes(self, **kwargs):
default is False
:type overwrite_energies: bool
:arg percentile: a percentile threshold to remove outliers, i.e. only showing data within *p*-th
to *100-p*-th percentile. Default is None, so no axis limits.
:type percentile: float
:arg vmin: a minimum value threshold to remove outliers, i.e. only showing data greater than vmin
This overrides percentile. Default is None, so no axis limits and little padding at the
bottom when energy=True.
:type vmin: float
:arg vmax: a maximum value threshold to remove outliers, i.e. only showing data less than vmax
This overrides percentile. Default is None, so no axis limits and a little padding for
interaction type labels.
:type vmax: float
'IB_solv' and 'IB_nosolv' are derived from empirical potentials from
O Keskin, I Bahar and colleagues from [OK98]_.
O Keskin, I Bahar and colleagues from [OK98]_ and have RT units.
'CS' is from MD simulations of amino acid pairs from Carlos Simmerling
and Gary Wu.
.. [OK98] Keskin O, Bahar I, Badretdinov AY, Ptitsyn OB, Jernigan RL,
Empirical solvent-mediated potentials hold for both intra-molecular
and inter-molecular inter-residue interactions
*Protein Sci* **1998** 7(12):2578-2586.
and Gary Wu for the InSty paper (under preparation) and have units kcal/mol.
"""

import numpy as np
Expand All @@ -3833,6 +3843,19 @@ def showCumulativeInteractionTypes(self, **kwargs):
if not isinstance(energy, bool):
raise TypeError('energy should be True or False')

p = kwargs.pop('percentile', None)
vmin = vmax = None
if p is not None:
vmin = np.percentile(matrix, p)
vmax = np.percentile(matrix, 100-p)

vmin = kwargs.pop('vmin', vmin)
vmax = kwargs.pop('vmax', vmax)

selstr = kwargs.pop('selstr', None)
if selstr is not None:
atoms = atoms.select(selstr)

ResNumb = atoms.select('protein and name CA').getResnums()
ResName = atoms.select('protein and name CA').getResnames()
ResChid = atoms.select('protein and name CA').getChids()
Expand All @@ -3856,16 +3879,27 @@ def showCumulativeInteractionTypes(self, **kwargs):
fig, ax = plt.subplots(num=None, figsize=(20,6), facecolor='w')
matplotlib.rcParams['font.size'] = '24'

if selstr is not None:
matrix_en_sum = sliceAtomicData(matrix_en_sum,
self._atoms.ca, selstr)

zeros_row = np.zeros(matrix_en_sum.shape)
pplot(zeros_row, atoms=atoms.ca)
pplot(zeros_row, atoms=atoms.ca, **kwargs)

ax.bar(ResList, matrix_en_sum, width, color='blue')

#plt.xlim([ResList[0]-0.5, ResList[-1]+0.5])
plt.ylim([min(matrix_en_sum)-1,0])
if vmin is None:
vmin = np.min(matrix_en_sum) * 1.2

if vmax is None:
vmax = np.max(matrix_en_sum)

plt.ylim([vmin, vmax])
plt.tight_layout()
plt.xlabel('Residue')
plt.ylabel('Cumulative Energy [kcal/mol]')

unit = 'RT' if energy_list_type in ['IB_solv', 'IB_nosolv'] else 'kcal/mol'
plt.ylabel('Cumulative Energy [{0}]'.format(unit))
plt.show()

return matrix_en_sum
Expand Down Expand Up @@ -3898,6 +3932,23 @@ def showCumulativeInteractionTypes(self, **kwargs):
matrix_hph_sum = np.sum(matrix_hph, axis=0)
matrix_dibs_sum = np.sum(matrix_dibs, axis=0)

all_ca = self._atoms.ca
if selstr is not None:
matrix_hbs_sum = sliceAtomicData(matrix_hbs_sum,
all_ca, selstr)
matrix_sbs_sum = sliceAtomicData(matrix_sbs_sum,
all_ca, selstr)
matrix_rib_sum = sliceAtomicData(matrix_rib_sum,
all_ca, selstr)
matrix_pistack_sum = sliceAtomicData(matrix_pistack_sum,
all_ca, selstr)
matrix_picat_sum = sliceAtomicData(matrix_picat_sum,
all_ca, selstr)
matrix_hph_sum = sliceAtomicData(matrix_hph_sum,
all_ca, selstr)
matrix_dibs_sum = sliceAtomicData(matrix_dibs_sum,
all_ca, selstr)

width = 0.8
fig, ax = plt.subplots(num=None, figsize=(20,6), facecolor='w')
matplotlib.rcParams['font.size'] = '24'
Expand Down Expand Up @@ -3940,7 +3991,14 @@ def showCumulativeInteractionTypes(self, **kwargs):
self._interactions_matrix = matrix_all

ax.legend(ncol=7, loc='upper center')
plt.ylim([0,max(sum_matrix)+3])

if vmin is None:
vmin = np.min(sum_matrix)

if vmax is None:
vmax = np.max(sum_matrix) * 1.5

plt.ylim([vmin, vmax])
plt.tight_layout()
plt.xlabel('Residue')
plt.ylabel('Number of counts')
Expand Down
8 changes: 8 additions & 0 deletions prody/utilities/catchall.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -610,6 +610,14 @@ def showMatrix(matrix, x_array=None, y_array=None, **kwargs):
to *100-p*-th percentile
:type percentile: float
:arg vmin: a minimum value threshold to remove outliers, i.e. only showing data greater than vmin
This overrides percentile.
:type vmin: float
:arg vmax: a maximum value threshold to remove outliers, i.e. only showing data less than vmax
This overrides percentile.
:type vmax: float
:arg interactive: turn on or off the interactive options
:type interactive: bool
Expand Down

0 comments on commit 3e3d822

Please sign in to comment.