import math
import os
import tempfile
import warnings
from typing import List
import matplotlib.gridspec as gridspec
import numpy
from matplotlib import artist
from matplotlib import figure
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from matplotlib.ticker import MaxNLocator
import schrodinger.application.desmond.report_helper as rh
import schrodinger.ui.qt.structure2d as structure2d
from schrodinger.application.desmond.constants import FEP_TYPES
from schrodinger.application.desmond.constants import ORIGINAL_INDEX
from schrodinger.infra import canvas2d
from schrodinger.Qt import QtGui
from schrodinger.ui import sketcher
from schrodinger.ui.qt.structure2d import CircleAnnotator
from schrodinger.ui.qt.structure2d import ColoredArrowAnnotator
from schrodinger.utils import fileutils
IGNORE_HYDROGEN = canvas2d.ChmAtomOption.H_Never
FIXUP = True
RESCALE = canvas2d.ChmCoordGenOption.Physical
STEREO = canvas2d.ChmMmctAdaptor.StereoFromGeometry_Safe
OPACITY = "opacity_value"
TORSIONS_PER_PAGE = 10
[docs]class InconsistentDataException(Exception):
pass
[docs]def create_sketcher():
sk = sketcher.SelectOnlySketcher()
sk.setLIDMode(True)
sk.setMinimizeWithoutIntersectingInteractions(False)
sk.menuBar().hide()
sk.LIDMainToolBar().hide()
sk.showBondCustomLabels(True)
# make interaction transparency depend on a property
style_keys = [
sketcher.sketcherStyle.H_BOND_TO_SIDE_CHAIN(),
sketcher.sketcherStyle.H_BOND_TO_BACKBONE(),
sketcher.sketcherStyle.METAL_INTERACTION(),
sketcher.sketcherStyle.PI_PI_INTERACTION(),
sketcher.sketcherStyle.PI_CAT_INTERACTION(),
sketcher.sketcherStyle.SALT_BRIDGE()
]
for style_key in style_keys:
style = sk.getScene().getStyle(style_key)
style.setFloatFunction("opacity", "MAP_FLOAT_PROPERTY", OPACITY + " 1")
return sk
[docs]def parse_res_tag(label):
try:
chainName, res_str, atomName = label.split(':')
except:
atomName = ''
chainName, res_str = label.split(':')
resName, resID = res_str.split('_')
return chainName, resName, resID, atomName
[docs]def get_residue_label(rname, resid, cname):
rl = rname + '\n' + resid
if cname != '_':
rl = cname + ':\n' + rl
return rl
[docs]def get_opacity(opacity_val):
o = opacity_val
if o > 1.0:
o = 1.0
elif o < 0.3:
o = 0.3
return o
def _new_sketcher_residue(sk, label, resid, resname, xyz):
res = sk.addResidue()
res.setLabel(label)
res.setNumber(resid)
res.setResidueType(resname)
res.set3DCoords(xyz[0], xyz[1], xyz[2])
return res
PL_INTERACTIONS = {
'hb': ('#7FC97F', 'H-bonds'),
'hphb': ('#BEAED4', 'Hydrophobic'),
'ion': ('#F0027F', 'Ionic'),
'wb': ('#386CB0', 'Water bridges')
}
ALL_PL_INTERACTIONS = ['hb', 'hphb', 'ion', 'wb']
[docs]class FEPEdgeReportMakerBase:
"""
Base class for all remport makers
"""
VERBOSE = False
[docs] def __init__(self, fep_edge_data, basename=None, perturbation_type=None):
"""
This a base class for generating various types of reports.
:type fep_edge_data: `FEPEdgeData`
:param fep_edge_data: Object containing all the data for this report
:type basename: string
:param basename: the basename of the file of the PDF report
:type perturbation_type: str
:param perturbation_type: FEP_TYPE value
"""
self.basename = basename
self.Elements = []
self.cleanup_tmp_file_list = []
self.schrod_tmp_dir = None
self.fed = fep_edge_data
self.perturbation_type = perturbation_type
def _load_modules(self):
self.app = rh.load_gui()
def _cleanup_temp_files(self):
if self.VERBOSE:
print("Cleaning up image directory:", self.schrod_tmp_dir)
if not self.schrod_tmp_dir:
return
fileutils.force_rmtree(self.schrod_tmp_dir, ignore_errors=True)
self.cleanup_tmp_file_list = []
self.schrod_tmp_dir = None
def _get_2d_ligand_image(self, st):
if not st:
return None
temp_file = self._get_temp_image_fn()
size = rh.generate_ligand_2d_image(temp_file, st, ret_size=True)
(sx, sy) = rh.aspectScale(size[0], size[1], 3., 2.25)
ligand_img = rh.platypus.Image(temp_file, sx * rh.inch, sy * rh.inch)
return ligand_img
[docs] def get_2d_ligand_image(self, lig_st):
try:
return self._get_2d_ligand_image(lig_st)
except:
return 'Cannot generate 2d'
def _get_temp_image_fn(self, prefix=''):
if not self.schrod_tmp_dir:
tmp_dir = fileutils.get_directory_path(fileutils.TEMP)
self.schrod_tmp_dir = tempfile.mkdtemp(dir=tmp_dir)
temp_fn = fileutils.get_next_filename(
os.path.join(self.schrod_tmp_dir, prefix + "image_tmp.png"), "")
self.cleanup_tmp_file_list.append(temp_fn)
return temp_fn
def _gen_free_energy_convergence_plot(self, start_time, end_time,
dG_forward, dG_forward_err,
dG_reverse, dG_reverse_err,
dG_sliding, dG_sliding_err):
fig = figure.Figure(figsize=(7.5, 2.0))
FigureCanvas(fig)
self._gen_free_energy_convergence_fig(start_time, end_time, dG_forward,
dG_forward_err, dG_reverse,
dG_reverse_err, dG_sliding,
dG_sliding_err, fig)
temp_plot = self._get_temp_image_fn()
with warnings.catch_warnings():
# DESMOND-10502
warnings.filterwarnings(
action="ignore",
category=UserWarning,
message=".*Axes that are not compatible with tight_layout",
)
fig.savefig(temp_plot, bbox_inches='tight', dpi=300)
size = fig.get_size_inches()
(sx, sy) = rh.aspectScale(size[0], size[1], 7.5, 2.5)
img = rh.platypus.Image(temp_plot, sx * rh.inch, (sy + .2) * rh.inch)
return img
def _get_free_energy_convergence_text(self, table=True):
text = '''
The total free energy differences between the two ligands (ΔG in
<i>kcal/mol</i>) in %s and %s legs are plotted as a function
of time. Three plots for each leg show the accumulated data during
different time window schemes: forward; reverse; and sliding window.
''' % (self.fed.leg1_name.lower(), self.fed.leg2_name.lower())
if table:
text += '''
<br/>
The tables report the associated bootstrap and analytical errors
estimates from corresponding simulation legs.
'''
return text
def _gen_free_energy_convergence_fig(self,
start_time,
end_time,
dG_forward,
dG_forward_err,
dG_reverse,
dG_reverse_err,
dG_sliding,
dG_sliding_err,
fig,
for_print=True):
global y_min
global y_max
y_min = 1000000
y_max = -1000000
def fill_error(ax, x, y, ye):
global y_min
global y_max
y_arr = numpy.array(y)
ye_arr = numpy.array(ye)
min_fill = y_arr - ye_arr
max_fill = y_arr + ye_arr
ax.fill_between(x,
y_arr - ye_arr,
y_arr + ye_arr,
facecolor='#888888',
edgecolor='None')
y_min = min(y_min, min(min_fill))
y_max = max(y_max, max(max_fill))
fig.clear()
fig.set_tight_layout(True)
time = numpy.linspace(start_time, end_time, len(dG_forward))
try:
timestep = time[1] - time[0]
except IndexError:
timestep = 0.3
ax_frw = fig.add_subplot(131)
ax_rev = fig.add_subplot(132, sharey=ax_frw)
ax_sld = fig.add_subplot(133, sharex=ax_frw, sharey=ax_frw)
ax_frw.plot(time, dG_forward, color='#990000')
ax_frw.set_xlabel("Simulation Time (nsec)", size='small')
ax_frw.set_ylabel(r"$\Delta$G (kcal/mol)", size='small')
fill_error(ax_frw, time, dG_forward, dG_forward_err)
# reverse plot should start with end time
dG_reverse.reverse()
dG_reverse_err.reverse()
time = numpy.linspace(start_time, end_time, len(dG_reverse))
ax_rev.plot(time, dG_reverse, color='#990000')
fill_error(ax_rev, time, dG_reverse, dG_reverse_err)
ax_rev.invert_xaxis()
ax_rev.set_xlabel("Reverse Simulation Time (nsec)", size='small')
artist.setp(ax_rev.get_yticklabels(), visible=False)
time = numpy.linspace(start_time, end_time, len(dG_sliding))
ax_sld.plot(time, dG_sliding, color='#990000')
fill_error(ax_sld, time, dG_sliding, dG_sliding_err)
ax_sld.set_xlabel("Sliding Time (nsec)", size='small')
ax_sld.set_ylabel(r"$\Delta$G (kcal/mol)",
size='small',
rotation=270,
labelpad=14)
ax_sld.yaxis.tick_right()
ax_sld.yaxis.set_label_position("right")
ax_frw.set_ylim((y_min, y_max))
ax_frw.set_xlim((start_time, end_time + timestep))
ax_rev.set_xlim((end_time + timestep, start_time))
if for_print:
for ax in [ax_frw, ax_rev, ax_sld]:
rh.change_plot_colors(ax)
[docs] def get_rest_density_img(self,
rest_density_data,
legend=True,
size_ratio=1.0):
ax_hist, fig = self.get_rest_density_plot(rest_density_data, legend,
True)
FigureCanvas(fig)
label_size = 8 if len(rest_density_data) <= 24 else 5
rh.change_plot_colors(ax_hist, label_size=label_size)
temp_plot = self._get_temp_image_fn()
fig.savefig(temp_plot, bbox_inches='tight', dpi=300)
size = fig.get_size_inches()
(sx, sy) = rh.aspectScale(size[0], size[1], 7.0 * size_ratio,
3.5 * size_ratio)
img = rh.platypus.Image(temp_plot, sx * rh.inch,
(sy + 0.5 * size_ratio) * rh.inch)
return img
[docs] def get_rest_density_plot(self,
rest_density_data,
legend=True,
for_print=False):
fig = figure.Figure(figsize=(5, 3.5))
ax_hist = fig.add_subplot(111)
nreplicas = len(rest_density_data)
bottom = numpy.zeros(nreplicas)
alpha = 0.8 if for_print else None
for irep, replica in enumerate(rest_density_data):
name = rh.replica_name(irep)
dist = numpy.array(replica)
color = rh.colors[irep % len(rh.colors)]
ax_hist.bar(list(range(nreplicas)),
dist,
label=name,
color=color,
width=1,
edgecolor=None,
linewidth=0,
align='center',
bottom=bottom,
alpha=alpha)
bottom += dist
ax_hist.set_xlim([0 - 0.5, nreplicas - 0.5])
ax_hist.set_xticks(list(range(nreplicas)))
ax_hist.set_xticklabels(list(range(nreplicas)),
fontsize=8 if for_print else None,
rotation=0 if nreplicas <= 24 else 70)
ax_hist.set_ylim([0, max(bottom)])
ax_hist.set_yticklabels([])
ax_hist.set_xlabel(r'$\lambda$ windows', size=8 if for_print else None)
if legend:
ax_hist.legend(loc=8,
bbox_to_anchor=(0.5, 1.02),
ncol=16 if nreplicas <= 24 else 20,
borderaxespad=0.,
borderpad=0.3,
handlelength=1,
handletextpad=0.2,
title='Replica Names',
columnspacing=0.4,
fancybox=True,
prop={'size': 7 if nreplicas <= 24 else 5})
return ax_hist, fig
@property
def get_rest_density_text(self):
text = """ \
For all the legs in the FEP simulation, each replica is color-coded and
the plot shows how it occupies different lambda windows during the
course of the simulation. Ideally, each replica will sample all lambda
windows uniformly, however, non-uniform sampling is sufficient in most
instances.
"""
return text
[docs] def get_torsions_plot(self,
tors1,
tors2,
tors_from,
tors_to,
for_print=False):
fig = figure.Figure()
FigureCanvas(fig)
# This is for generating PDF report, no interaction needed.
self.create_torsions_plot(fig,
tors1,
tors2,
tors_from,
tors_to,
None,
None,
None,
for_print=for_print)
return fig
def _align_chmmol(self,
ligand_chmmol,
ligand_core_idx=None,
template_chmmol=None,
template_core_idx=None):
decrement = lambda x: x - 1
if template_chmmol is not None:
ligand_core_idx = list(map(decrement, ligand_core_idx))
template_core_idx = list(map(decrement, template_core_idx))
canvas2d.Chm2DCoordGen.generateFromTemplateAndApply(
ligand_chmmol, template_chmmol, ligand_core_idx,
template_core_idx, IGNORE_HYDROGEN, RESCALE, FIXUP)
else:
canvas2d.Chm2DCoordGen.generateAndApply(
ligand_chmmol, IGNORE_HYDROGEN,
canvas2d.ChmCoordGenOption.Rendering)
[docs] def get_structure_item(self, st, rect=None):
"""
Returns structure as a graphic item to draw and to annotate.
:param st: Structure to draw in 2D
:type st: `structure.Structure`
:param rect: Rectangle (bounding box) for the structure
:type rect: QtCore.QRectF
:return: Structure item to be drawn in 2D.
:rtype: `structure2d.structure_item`
"""
si = structure2d.structure_item(rect)
si.model2d.setShowHydrogenPreference(0)
si.set_structure(st)
return si
[docs] def get_2d_tors_annotated_lig_pair(self,
tors1,
tors2,
tors_from,
tors_to,
rect=None):
"""
Returns structure scenes for two 2d ligands with annotations.
"""
def process_tor(itor, tor, color):
if itor < len(tor):
a = tor[itor]._sol_atoms[1:3]
return [a[0], a[1], False, color]
def add_annotators(structure_item, bondlist):
structure_item.clear_annotators()
annotator = ColoredArrowAnnotator(structure_item.chmmol,
bond_info=bondlist,
draw_arrowhead=False)
structure_item.add_annotator(annotator)
structure_item.model2d.setColorMap(rh.blackColorMap)
structure_item.generate_picture(gen_coord=False)
bondlist1 = []
bondlist2 = []
for itor in range(tors_from, tors_to):
color = rh.get_qcolor(rh.colors[itor % len(rh.colors)])
b1 = process_tor(itor, tors1, color)
b2 = process_tor(itor, tors2, color)
if b1 is not None:
bondlist1.append(b1)
if b2 is not None:
bondlist2.append(b2)
st1, st2 = self._getTorsionStructures()
lig1_item = self.get_structure_item(st1, rect)
lig2_item = self.get_structure_item(st2, rect)
self._align_chmmol(lig1_item.chmmol)
self._align_chmmol(lig2_item.chmmol, self.fed._atom_mapping[1],
lig1_item.chmmol, self.fed._atom_mapping[0])
add_annotators(lig1_item, bondlist1)
add_annotators(lig2_item, bondlist2)
return lig1_item, lig2_item
def _get_annotated_2d_lig_img(self, lig_2d_si):
temp_img_fn = self._get_temp_image_fn()
try:
size = rh.save_2d_annotated_img(lig_2d_si,
temp_img_fn,
crop=True,
ret_size=True)
(sx, sy) = rh.aspectScale(size[0], size[1], 3.25, 2.25)
img = rh.platypus.Image(temp_img_fn, sx * rh.inch, sy * rh.inch)
except:
img = '2D image cannot be generated'
return img
[docs]class FEPEdgeReportMaker(FEPEdgeReportMakerBase):
HELIX = 'helix'
STRAND = 'strand'
HELIX_COLOR = '#F06040'
STRAND_COLOR = '#40C0E0'
B_FACTOR = '#9400D3'
[docs] def __init__(self,
fep_edge_data,
basename=None,
perturbation_type=FEP_TYPES.SMALL_MOLECULE):
"""
This class generates a PDF report for an FEP/REST (+) type of job. Both
Legs of the simulations are processed in the `FEPEdgeData` and is used
by this class to compile a result.
:type fep_edge_data: `FEPEdgeData`
:param fep_edge_data: Object containing all the data for this report
:type basename: string
:param basename: the basename of the file of the PDF report
:type perturbation_type: str
:param perturbation_type: Type of FEP job it was
"""
super(FEPEdgeReportMaker, self).__init__(fep_edge_data, basename,
perturbation_type)
self._load_modules()
# The following are needed for handling mouse hovering on torsions plot.
# In the following dictionaries, the key is page index of torsions.
self.fig_dict = {}
self.gs_dict = {}
self.lig1_item_dict = {}
self.lig2_item_dict = {}
self.tors1_dict = {}
self.tors2_dict = {}
self.torsion_colors_dict = {}
self.torsions_page_index = 0
def _load_modules(self):
self.app = rh.load_gui()
from schrodinger.ui.sequencealignment.sequence_viewer import \
SequenceViewer
self.SeqViewer = SequenceViewer
from schrodinger.ui.sequencealignment.globals import ANNOTATION_RESNUM
from schrodinger.ui.sequencealignment.globals import ANNOTATION_SSBOND
from schrodinger.ui.sequencealignment.globals import COLOR_POSITION
self.COLOR_POSITION = COLOR_POSITION
self.ANNOTATION_SSBOND = ANNOTATION_SSBOND
self.ANNOTATION_RESNUM = ANNOTATION_RESNUM
[docs] def report(self):
self._init_report()
self._report_fep_simulation_details()
self._ligand_perturbation_details()
self._protein_details()
self._salt_details()
self._free_energy_convergence_profiles()
self._replica_exchange_density()
self._protein_report()
self._ligand_report()
self._ligand_torsion_report()
self._protein_ligand_report()
self.doc.build(self.Elements, canvasmaker=rh.NumberedCanvas)
self._cleanup_temp_files()
print(self.basename, 'report is written')
def _init_report(self):
pdf_filename = self.basename + '.pdf'
self.doc = rh.platypus.SimpleDocTemplate(pdf_filename)
self.doc.title = u'FEP+ Report by Schrodinger Inc.'
self.doc.author = "Dmitry Lupyan, PhD"
self.doc.leftMargin = 30.
self.doc.rightMargin = 20.
self.doc.topMargin = 10.
self.doc.bottomMargin = 10.
rh.report_add_logo(self.Elements)
rh.header(self.Elements, "FEP+ Simulation Details and Results")
def _replica_exchange_density(self):
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Exchange Density of FEP Replicas Over "
"λ-Windows")
sol_exchanges = fed.sol_rest_exchanges
cpx_exchanges = fed.cpx_rest_exchanges
rh.pargph(Ele, '<u>%s Leg</u>' % fed.leg1_name)
rh.add_spacer(Ele)
sol_img = self.get_rest_density_img(sol_exchanges)
Ele.append(sol_img)
rh.pargph(Ele, '<u>%s Leg</u>' % fed.leg2_name)
rh.add_spacer(Ele)
cpx_img = self.get_rest_density_img(cpx_exchanges)
Ele.append(cpx_img)
rh.add_spacer(Ele)
rh.pargph(Ele, self.get_rest_density_text, fontsize=10)
@property
def get_rest_density_text(self):
text = """ \
For all the legs in the FEP simulation, each replica is color-coded and
the plot shows how it occupies different lambda windows during the
course of the simulation. Ideally, each replica will sample all lambda
windows uniformly, however, non-uniform sampling is sufficient in most
instances.
"""
return text
def _protein_ligand_report(self):
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(
Ele, "Protein-Ligand Interactions for End-Point "
"λ-Replicas")
residues = fed.cpx_sid_protein_residues
pli0 = fed.cpx_sid_pl_results[0]
pli1 = fed.cpx_sid_pl_results[1]
nframes = fed.cpx_sid_number_of_frames
bar_chart = self._get_plot_pl_bar_chart(residues, pli0, pli1, nframes)
Ele.append(bar_chart)
rh.pargph(Ele, self._get_pl_text(), fontsize=10)
lp_stats = fed.cpx_sid_lp_results
cutoff = 0.20
lig_noh_atom_dict, lig_atom_dict = fed.get_ligand1_atom_dict()
atom_mapping = fed.atom_mapping
lig1_st, lig2_st = fed.ligand1_st, fed.ligand2_st
if fed.perturbation_type == FEP_TYPES.COVALENT_LIGAND:
lig1_st, lig2_st = fed.ligand1_cpx_st, fed.ligand2_cpx_st
lid0_img = self._get_lid_img(lp_stats[0],
lig_noh_atom_dict,
lig_atom_dict,
lig1_st,
atom_mapping[0],
template_st=None,
template_core_idx=None,
cutoff=cutoff)
lig_noh_atom_dict, lig_atom_dict = fed.get_ligand2_atom_dict()
lid1_img = self._get_lid_img(lp_stats[1],
lig_noh_atom_dict,
lig_atom_dict,
lig2_st,
atom_mapping[1],
template_st=None,
template_core_idx=None,
cutoff=cutoff)
table = [['Ligand 1', 'Ligand 2'], [lid0_img, lid1_img]]
width = [4.0, 4.0]
nfields = len(table[0])
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
def _get_lid_img(self,
stats,
lig_noh_dict,
lig_atom_dict,
ligand_st,
ligand_core_idx,
template_st=None,
template_core_idx=None,
cutoff=0.2):
lid = self._get_lid(stats, lig_noh_dict, lig_atom_dict, ligand_st,
ligand_core_idx, template_st, template_core_idx,
cutoff)
lid_image = lid.getQImage(2400, 2400)
del lid
image_fn = self._get_temp_image_fn()
lid_image = rh.crop_image(lid_image)
lid_image.save(image_fn)
size = lid_image.size
(sx, sy) = rh.aspectScale(size[0], size[1], 3.75, 4.0)
return rh.platypus.Image(image_fn, sx * rh.inch, sy * rh.inch)
@staticmethod
def _st_to_chmmol(ligand_st):
"""
Convert ligand_st to Canvas 2d structure
:param ligand_st: Ligand structure
:type ligand_st: `schrodinger.structure`
:returns: Canvas molecule
:rtype: `canvas2d.ChmMol`
"""
adaptor = canvas2d.ChmMmctAdaptor()
return adaptor.create(ligand_st.handle, STEREO)
def _get_aligned_chmmol(self,
ligand_st,
ligand_core_idx,
template_st=None,
template_core_idx=None):
ligand_chmmol = self._st_to_chmmol(ligand_st)
if template_st is None:
template_chmmol = ligand_chmmol
template_core_idx = ligand_core_idx
else:
template_chmmol = self._st_to_chmmol(template_st)
self._align_chmmol(ligand_chmmol, ligand_core_idx, template_chmmol,
template_core_idx)
return ligand_chmmol
def _get_lid(self,
stats,
lig_noh_dict,
lig_atom_dict,
ligand_st,
ligand_core_idx,
template_st=None,
template_core_idx=None,
cutoff=0.2):
sk = create_sketcher()
ligand_chmmol = self._get_aligned_chmmol(ligand_st, ligand_core_idx,
template_st, template_core_idx)
center = False
regenerateCoords = False
ligand_atoms = sk.addLigand(ligand_chmmol, center, regenerateCoords)
#copy over 3d coordinates (only heavy atoms are kept in the 2d molecule)
for i, a in enumerate(
[a for a in ligand_st.atom if a.atomic_number != 1]):
ligand_atoms[i].set3DCoords(*a.xyz)
already_added_res_dict = {}
already_added_metal_dict = {}
for hbond in stats['hbond']:
if hbond[0] > cutoff:
self._lid_add_hbonds(hbond, sk, already_added_res_dict,
ligand_st, lig_atom_dict, lig_noh_dict,
ligand_atoms)
for ionic in stats['ionic']:
if ionic[0] > cutoff:
self._lid_add_ionic(ionic, sk, already_added_res_dict,
lig_noh_dict, ligand_atoms)
for intra_hb in stats['l_hbond']:
if intra_hb[0] > cutoff:
self._lid_add_intramolecular_hb(intra_hb, sk, ligand_st,
lig_atom_dict, lig_noh_dict,
ligand_atoms)
for l_metal in stats['l_metal']:
if l_metal[0] > cutoff:
self._lid_add_lig_metal(l_metal, sk, already_added_metal_dict,
lig_noh_dict, ligand_atoms)
for p_metal in stats['p_metal']:
if p_metal[0] > cutoff:
self._lid_add_prot_metal(p_metal, sk, already_added_metal_dict,
already_added_res_dict)
for pipi in stats['pipi']:
if pipi[0] > cutoff:
self._lid_add_pipi(pipi, sk, already_added_res_dict,
lig_noh_dict, ligand_atoms)
for picat in stats['picat']:
if picat[0] > cutoff:
self._lid_add_picat(picat, sk, already_added_res_dict,
lig_noh_dict, ligand_atoms)
for hphob in stats['hphob']:
if hphob[0] > cutoff:
self._lid_add_hydrophobic(hphob, sk, already_added_res_dict)
for wat_br in stats['wat_br']:
if wat_br[0] > cutoff:
self._lid_add_water_br(wat_br, sk, already_added_res_dict,
ligand_st, lig_atom_dict, lig_noh_dict,
ligand_atoms)
for lw in stats['lig_wat']:
if lw[0] > 0:
self._lid_add_lig_water_exposed(lw, sk, ligand_st, lig_noh_dict,
ligand_atoms)
sk.cleanUp(True)
return sk
@staticmethod
def _lid_add_lig_water_exposed(res_label, sk, ligand_st, lig_noh_dict,
ligand_atoms):
val = res_label[0] * 3
if val > 15.0:
val = 15.
for aid in res_label[1]:
if aid in lig_noh_dict:
if val > 5.0:
index = lig_noh_dict[aid]
sk.setSolventExposureForAtom(index - 1, val)
@staticmethod
def _lid_add_water_br(res_label, sk, already_added_res_dict, ligand_st,
lig_atom_dict, lig_noh_dict, ligand_atoms):
"""
`res_label` for water-mediated interactions looks like this::
(0.610, 'B:LEU_52:O', [-7.148, -19.094, -5.2140861], 3702)
(freq, prot_label, prot_xyz, ligand_aid)
"""
cname, rname, resid, atom_name = parse_res_tag(res_label[1])
rl = get_residue_label(rname, resid, cname)
if rl not in already_added_res_dict:
res = _new_sketcher_residue(sk, rl, int(resid), rname, res_label[2])
already_added_res_dict[rl] = res
else:
res = already_added_res_dict[rl]
lig_atm_idx = lig_atom_dict[res_label[3]]
lig_ele = ligand_st.atom[lig_atm_idx].element
water = sk.addResidue()
water.setLabel('H2O')
water.setResidueType('H2O')
if lig_ele == 'H':
# get index of the heavy atom that's bonded to the hydrogen
bond = ligand_st.atom[lig_atm_idx].bond[1]
lig_atm_idx = lig_noh_dict[bond.atom2.property[ORIGINAL_INDEX]]
hb = sk.addResidueInteraction(ligand_atoms[lig_atm_idx - 1], water)
else:
lig_atm_idx = lig_noh_dict[res_label[3]]
hb = sk.addResidueInteraction(water, ligand_atoms[lig_atm_idx - 1])
sk.setSolventExposureForAtom(lig_atm_idx - 1, 5.01)
if atom_name[0:1] == 'O':
hb2 = sk.addResidueInteraction(water, res)
else:
hb2 = sk.addResidueInteraction(res, water)
style_key = sketcher.sketcherStyle.H_BOND_TO_SIDE_CHAIN()
backbone_atoms = ['H', 'O', 'HA2', 'HA3']
if atom_name in backbone_atoms:
style_key = sketcher.sketcherStyle.H_BOND_TO_BACKBONE()
hb.setStyleKey(sketcher.sketcherStyle.H_BOND_TO_BACKBONE())
hb.setFloatProperty(OPACITY, get_opacity(res_label[0]))
hb2.setStyleKey(style_key)
hb2.setFloatProperty(OPACITY, get_opacity(res_label[0]))
label = str(int(res_label[0] * 100)) + '%'
hb2.setCustomLabel(label)
hb.setCustomLabel(label)
@staticmethod
def _lid_add_hydrophobic(res_label, sk, already_added_res_dict):
"""
`res_label` for non-specific hydrophobic interactions looks like this::
(0.801, 'B:LEU_52', [-5.9, -21.1, -5.4],
[3694, 3695, 3696, 3709, 3710, 3711])
(freq, prot_label, [prot_xyz], fragment_aids)
"""
cname, rname, resid, atom_name = parse_res_tag(res_label[1])
rl = get_residue_label(rname, resid, cname)
if rl not in already_added_res_dict:
res = _new_sketcher_residue(sk, rl, int(resid), rname, res_label[2])
already_added_res_dict[rl] = res
@staticmethod
def _lid_add_picat(res_label, sk, already_added_res_dict, lig_atom_dict,
ligand_atoms):
"""
`res_label` for pi-cat interactions looks like this::
(0.16, 'B:ARG_77:CZ', [-5.9, -11.1, -5.9], [698, 699, 701, 702, 700])
(freq, prot_label, [prot_xyz], lig_ring_aids)
in case the aromatic group is on the protein, and charged group on the
ligand::
(1.0, 'B:TYR_99', [13.404797, 3.522047, 15.211519], 2945)
(freq, prot_label (aromatic) , [prot_xyz], ligand_aid)
"""
cname, rname, resid, atom_name = parse_res_tag(res_label[1])
rl = get_residue_label(rname, resid, cname)
if rl not in already_added_res_dict:
res = _new_sketcher_residue(sk, rl, int(resid), rname, res_label[2])
already_added_res_dict[rl] = res
else:
res = already_added_res_dict[rl]
if isinstance(res_label[3], list):
lig_ring_atm_idx = [lig_atom_dict[r] - 1 for r in res_label[3]]
picat = sk.addResidueInteraction(res,
ligand_atoms[lig_ring_atm_idx[0]])
atoms = [ligand_atoms[r] for r in lig_ring_atm_idx]
picat.setOtherEndAtoms(atoms)
else:
lig_atm_idx = lig_atom_dict[res_label[3]]
picat = sk.addResidueInteraction(res, ligand_atoms[lig_atm_idx - 1])
picat.setStyleKey(sketcher.sketcherStyle.PI_CAT_INTERACTION())
picat.setFloatProperty(OPACITY, get_opacity(res_label[0]))
label = str(int(res_label[0] * 100)) + '%'
picat.setCustomLabel(label)
@staticmethod
def _lid_add_pipi(res_label, sk, already_added_res_dict, lig_atom_dict,
ligand_atoms):
"""
`res_label` for pi-pi interactions looks like this::
(0.068, 'B:TYR_14', [4.0, -22.1, -8.7], [3698, 3699, 3701, 3702, 3700],
0.0)
(freq, prot_label, [Calpha_xyz], ligand_rings_aid, f2f_fraction)
"""
cname, rname, resid, atom_name = parse_res_tag(res_label[1])
rl = get_residue_label(rname, resid, cname)
if rl not in already_added_res_dict:
res = _new_sketcher_residue(sk, rl, int(resid), rname, res_label[2])
already_added_res_dict[rl] = res
else:
res = already_added_res_dict[rl]
lig_ring_atm_idx = []
for r in res_label[3]:
index = lig_atom_dict[r] - 1
lig_ring_atm_idx.append(index)
pipi = sk.addResidueInteraction(res, ligand_atoms[lig_ring_atm_idx[0]])
atoms = [ligand_atoms[r] for r in lig_ring_atm_idx]
pipi.setOtherEndAtoms(atoms)
if res_label[4] >= 0.5:
pipi.setIsFaceToFace(True)
pipi.setStyleKey(sketcher.sketcherStyle.PI_PI_INTERACTION())
pipi.setFloatProperty(OPACITY, get_opacity(res_label[0]))
label = str(int(res_label[0] * 100)) + '%'
pipi.setCustomLabel(label)
@staticmethod
def _lid_add_prot_metal(res_label, sk, already_added_metal_dict,
already_added_res_dict):
"""
`res_label` for metal-protein mediated interactions looks like this::
(1.0, 'A:ASP_116:OD1', 1511, [9.3, 17.1, -3.8], 1, [7.4, 16.6,
-5.2])
(freq, prot_label, Calpha_aid, [Calpha_xyz], metal_aid,
metal_label, metal_xyz)
"""
cname, rname, resid, atom_name = parse_res_tag(res_label[1])
rl = get_residue_label(rname, resid, cname)
if rl not in already_added_res_dict:
res = _new_sketcher_residue(sk, rl, int(resid), rname, res_label[3])
already_added_res_dict[rl] = res
else:
res = already_added_res_dict[rl]
# check if metal already has been added
metal_aid = res_label[4]
if metal_aid not in already_added_metal_dict:
mres = sk.addResidue()
lab = res_label[5].split('_')[1].split(':')[0]
xyz = res_label[6]
mres.setLabel(lab)
mres.setResidueType('Metal')
mres.set3DCoords(xyz[0], xyz[1], xyz[2])
already_added_metal_dict[metal_aid] = mres
else:
mres = already_added_metal_dict[metal_aid]
metal = sk.addResidueInteraction(res, mres)
metal.setStyleKey(sketcher.sketcherStyle.METAL_INTERACTION())
metal.setFloatProperty(OPACITY, get_opacity(res_label[0]))
label = str(int(res_label[0] * 100)) + '%'
metal.setCustomLabel(label)
@staticmethod
def _lid_add_lig_metal(res_label, sk, already_added_metal_dict,
lig_atom_dict, ligand_atoms):
"""
`res_label` for metal-ligand mediated interactions looks like this::
(1.0, 'L-FRAG_0:O1754', 1754, [9.4, 16.6, -7.0], 1, [7.4, 16.6,
-5.2])
(freq, lig_label, [ligand_xyz], metal_aid, metal_label, metal_xyz)
"""
metal_aid = res_label[4]
if metal_aid not in already_added_metal_dict:
mres = sk.addResidue()
lab = res_label[5].split('_')[1].split(':')[0]
xyz = res_label[6]
mres.setLabel(lab)
mres.setResidueType('Metal')
mres.set3DCoords(xyz[0], xyz[1], xyz[2])
already_added_metal_dict[metal_aid] = mres
else:
mres = already_added_metal_dict[metal_aid]
lig_atm_idx = lig_atom_dict[res_label[2]]
metal = sk.addResidueInteraction(ligand_atoms[lig_atm_idx - 1], mres)
metal.setStyleKey(sketcher.sketcherStyle.METAL_INTERACTION())
metal.setFloatProperty(OPACITY, get_opacity(res_label[0]))
label = str(int(res_label[0] * 100)) + '%'
metal.setCustomLabel(label)
@staticmethod
def _lid_add_intramolecular_hb(res_label, sk, ligand_st, lig_atom_dict,
lig_noh_dict, ligand_atoms):
"""
res_label for internal (intramolecular) hydrogen bonds looks like this:
(0.123, 4720, 4764)
freq, lig_donor_aid, lig_acceptor_aid
"""
donor_aid = lig_atom_dict[res_label[1]]
# get index of the heavy atom that's bonded to the hydrogen
heavy_aid = ligand_st.atom[donor_aid].bond[1].atom2.index
# get original index of the heavy atoms connected the hydrogen
donor_heavy_aid = ligand_st.atom[heavy_aid].property[ORIGINAL_INDEX]
donor_aid = lig_noh_dict[donor_heavy_aid]
accept_aid = lig_noh_dict[res_label[2]]
# add LID hbond
intra_hb = sk.addResidueInteraction(ligand_atoms[donor_aid - 1],
ligand_atoms[accept_aid - 1])
intra_hb.setStyleKey(sketcher.sketcherStyle.H_BOND_TO_SIDE_CHAIN())
intra_hb.setFloatProperty(OPACITY, get_opacity(res_label[0]))
label = str(int(res_label[0] * 100)) + '%'
intra_hb.setCustomLabel(label)
@staticmethod
def _lid_add_ionic(res_label, sk, already_added_res_dict, lig_noh_dict,
ligand_atoms):
"""
res_label for ionic bonds looks like this::
(0.0148, 'B:ARG_13:NH2', 180, [4.014, -2.630, 8.924], 2357)
freq, prot_label, prot_aid, [prot_ xyz], lig_aid
"""
cname, rname, resid, atom_name = parse_res_tag(res_label[1])
rl = get_residue_label(rname, resid, cname)
if rl not in already_added_res_dict:
res = _new_sketcher_residue(sk, rl, int(resid), rname, res_label[3])
already_added_res_dict[rl] = res
else:
res = already_added_res_dict[rl]
if res_label[4] in lig_noh_dict:
lig_atom_idx = lig_noh_dict.get(res_label[4])
else:
return
ionic = sk.addResidueInteraction(ligand_atoms[lig_atom_idx - 1], res)
ionic.setStyleKey(sketcher.sketcherStyle.SALT_BRIDGE())
ionic.setFloatProperty(OPACITY, get_opacity(res_label[0]))
label = str(int(res_label[0] * 100)) + '%'
ionic.setCustomLabel(label)
@staticmethod
def _lid_add_hbonds(res_label, sk, already_added_res_dict, ligand_st,
lig_atom_dict, lig_noh_dict, ligand_atoms):
"""
res_label for hydrogen bonds looks like this:
(0.40, 'A:LYS_97:HZ1', 2005, [-9.87, 7.35, -35.22], 'L-FRAG_0:O16', 16)
frequency, prot_label prot_aid, [xyz], lig_label, lig_aid
"""
cname, rname, resid, atom_name = parse_res_tag(res_label[1])
rl = get_residue_label(rname, resid, cname)
if rl not in already_added_res_dict:
res = _new_sketcher_residue(sk, rl, int(resid), rname, res_label[3])
already_added_res_dict[rl] = res
else:
res = already_added_res_dict[rl]
lig_atm_idx = lig_atom_dict[res_label[5]]
if ligand_st.atom[lig_atm_idx].element == 'H':
# get index of the atom that's bonded to the hydrogen
bond = ligand_st.atom[lig_atm_idx].bond[1]
lig_atm_idx = lig_noh_dict[
bond.atom2.property['i_m_original_index']]
hb = sk.addResidueInteraction(ligand_atoms[lig_atm_idx - 1], res)
else:
lig_atm_idx = lig_noh_dict[res_label[5]]
hb = sk.addResidueInteraction(res, ligand_atoms[lig_atm_idx - 1])
style_key = sketcher.sketcherStyle.H_BOND_TO_SIDE_CHAIN()
backbone_atoms = [
'H', 'O', 'HA', 'HA2', 'HA3', '1H', '2H', '3H', '1HA', '2HA', '3HA'
]
if atom_name in backbone_atoms:
style_key = sketcher.sketcherStyle.H_BOND_TO_BACKBONE()
hb.setStyleKey(style_key)
hb.setFloatProperty(OPACITY, get_opacity(res_label[0]))
label = str(int(res_label[0] * 100)) + '%'
hb.setCustomLabel(label)
def _get_pl_text(self):
text = """\
Above bar charts illustrate the type of interactions the protein
residues make with each ligand. Multiple types of specific
interactions are monitored throughout the simulation and provide a great
way to examine and compare how the protein interacts with both ligands.
These differences may account for the binding energy change between the
two ligands.
The specific interactions types monitored and displayed are:
<b><font color='#7FC97F'> hydrogen bond</font></b>,
<b><font color='#BEAED4'>hydrophobic</font></b>,
<b><font color='#F0027F'>ionic</font></b> and
<b><font color='#386CB0'>water bridges</font></b>. The stacked bar charts
are normalized over the course of the trajectory. More information
about the geometry and the different interaction subtype categories
can be found in the Schrödinger's <i>Desmond User Manual</i>,
under 'Simulation Interactions Diagram' (SID) section. <i><b>Note:</b></i>
The values may exceed 100% as the residue can make multiple contacts of
the same type with the ligand.
"""
return text
def _get_plot_pl_bar_chart(self,
label,
data_dict0,
data_dict1,
nframes,
name_1="Ligand 1",
name_2="Ligand 2"):
fig = self._create_pl_bar_fig(label,
data_dict0,
data_dict1,
nframes,
for_print=True,
name_1=name_1,
name_2=name_2)
FigureCanvas(fig)
temp_plot = self._get_temp_image_fn()
fig.savefig(temp_plot, bbox_inches='tight', dpi=300)
size = fig.get_size_inches()
(sx, sy) = rh.aspectScale(size[0], size[1], 6.0, 4.0)
img = rh.platypus.Image(temp_plot, sx * rh.inch, sy * rh.inch)
return img
def _create_pl_bar_fig(self,
label,
data_dict0,
data_dict1,
nframes,
pl_interactions=ALL_PL_INTERACTIONS,
name_1="Ligand 1",
name_2="Ligand 2",
fig=None,
for_print=False):
def check_if_one_chain(labels_list):
"""
Check if residue labels have the same chain name
"""
if len(set(lab[0] for lab in labels_list)) == 1:
return True
return False
# copy and reverse the labels, such that the first residue is at the top
label = label[::-1]
nresidues = len(label)
rrange = numpy.arange(nresidues)
sum_list0 = numpy.zeros(nresidues)
sum_list1 = numpy.zeros(nresidues)
gs = gridspec.GridSpec(1, 2, top=0.95, bottom=0.07, wspace=0.35)
if fig is None:
fig = figure.Figure(figsize=(6.0, 4.0))
else:
fig.clear()
ax1 = fig.add_subplot(gs[0])
ax2 = fig.add_subplot(gs[1])
for type in pl_interactions:
itype_vals0 = numpy.array([data_dict0[lab][type] for lab in label
]) / float(nframes) * 100
ax1.barh(rrange,
itype_vals0,
left=sum_list0,
color=PL_INTERACTIONS[type][0],
label=type)
sum_list0 += itype_vals0
itype_vals1 = numpy.array([data_dict1[lab][type] for lab in label
]) / float(nframes) * 100
ax2.barh(rrange,
itype_vals1,
left=sum_list1,
color=PL_INTERACTIONS[type][0],
label=type)
sum_list1 += itype_vals1
fontsize = 8 if for_print else 11
ax1.yaxis.tick_right()
ax1.yaxis.set_ticks(rrange + 0.5)
if check_if_one_chain(label):
new_label = [lab[2:].replace('_', ' ') for lab in label]
else:
new_label = [lab.replace('_', ' ') for lab in label]
ax1.yaxis.set_ticklabels(new_label, fontsize=fontsize)
ax2.yaxis.tick_left()
ax2.yaxis.set_ticks(rrange + 0.5)
xlim1 = ax1.get_xlim()
xlim2 = ax2.get_xlim()
# set same ranges for both plots
xlim_max = max(xlim1[1], xlim2[1]) + 5
ax1.set_xlim((0, xlim_max))
ax2.set_xlim((0, xlim_max))
ax1.set_ylim(-0.5, nresidues + 0.25)
ax2.set_ylim(-0.5, nresidues + 0.25)
ax1.invert_xaxis()
ax2.yaxis.set_ticklabels(len(label) * [], fontsize=fontsize)
xlabel_1 = '%% interact. w/ %s' % name_1
xlabel_2 = '%% interact. w/ %s' % name_2
ax1.set_xlabel(xlabel_1, fontsize=fontsize)
ax2.set_xlabel(xlabel_2, fontsize=fontsize)
if for_print:
for tick in ax1.get_xticklabels() + ax2.get_xticklabels():
tick.set_fontsize(7.5)
rh.change_plot_colors(ax1, ticks=False, labels=False)
rh.change_plot_colors(ax2, ticks=False, labels=False)
return fig
def _getTorsionStructures(self):
"""
:return: the two structures for the torsion plots
"""
return self.fed.ligand1_st, self.fed.ligand2_st
[docs] def remove_last_atom_highlight(self, lig_item):
"""
Removes previously highlighted atoms, if any, in a torsion.
:param lig_item: Ligand structure item.
:type lig_item: `structure2d.structure_item`
"""
if (lig_item.annotators and
isinstance(lig_item.annotators[-1], CircleAnnotator)):
del lig_item.annotators[-1]
lig_item.generate_picture(gen_coord=False)
[docs] def add_atom_highlight(self, lig_item, torsion, color):
"""
Add new highlight for atoms in the given torsion.
:param lig_item: Ligand structure item.
:type lig_item: `structure2d.structure_item`
:param torsion: Ligand torsion to highlight
:type torsion: `fep_edge_data.FEPTorsions`
:param color: Highlight color
:type color: `QtGui.QColor`
"""
self.remove_last_atom_highlight(lig_item)
anno = CircleAnnotator(torsion._sol_atoms[0],
color=color,
gradient=True)
anno.addAtom(torsion._sol_atoms[3])
lig_item.add_annotator(anno)
lig_item.generate_picture(gen_coord=False)
[docs] def on_move(self, event):
"""
The `mpl_connect` motion_notify_event slot for the torsion plot figure
canvas.
If the mouse is inside a torsion plot, highlight the 1st and 4th
atoms in that torsion in the structure.
Because torsions can spread to multiple tabs (or pages), need to
access figure and structure items for each tab.
:param event: The mouse event that triggered this callback
:type event: `matplotlib.backend_bases.MouseEvent`
"""
page = self.torsions_page_index
lig1_item = self.lig1_item_dict[page]
lig2_item = self.lig2_item_dict[page]
if event.xdata is None or event.ydata is None:
# outside the subplots
self.remove_last_atom_highlight(lig1_item)
self.remove_last_atom_highlight(lig2_item)
return
fig = self.fig_dict[page]
fig_size = fig.get_size_inches() * fig.dpi
fig_height = fig_size[1]
# grid positions are the bottom, top, left and right locations
# of each grid relative to the entire figure.
grid_positions = self.gs_dict[page].get_grid_positions(fig)
subplot_bottoms = grid_positions[0]
subplot_tops = grid_positions[1]
tors1 = self.tors1_dict[page]
tors2 = self.tors2_dict[page]
# The number of rows per page
nrow = len(subplot_bottoms)
# If the mouse is inside a torsion plot, atoms from the corresponding
# torsions from both ligands, if exist, should be highlighted. When
# the mouse is moved to a different torsion plot, highlight from previous
# atoms should be removed.
for r in range(nrow):
if (event.y >= subplot_bottoms[r] * fig_height and
event.y <= subplot_tops[r] * fig_height):
color = self.torsion_colors_dict[page][r]
self._highlight_matching_torsion_atoms(r, tors1, lig1_item,
color)
self._highlight_matching_torsion_atoms(r, tors2, lig2_item,
color)
return
def _highlight_matching_torsion_atoms(self, row, tors, lig_item, color):
"""
Highlight atoms in a torsion that matches the plot where
the mouse is in.
:param row: Row index of the torsion plot that the mouse is in
:type row: int
:param tors: List of torsions from a ligand
:type tors: list
:param lig_item: Structure item for a ligand.
:type lig_item: `structure2d.structure_item`
:param color: Color used to highlight the atoms.
:type color: `QtGui.QColor`
"""
if row < len(tors) and tors[row]:
self.add_atom_highlight(lig_item, tors[row], color)
else:
self.remove_last_atom_highlight(lig_item)
[docs] def create_torsions_plot(self,
fig,
tors1,
tors2,
tors_from,
tors_to,
ipage=None,
structure1_item=None,
structure2_item=None,
for_print=False):
"""
Creates a plot for torsions using the given matplot figure.
:param fig: Figure to draw the torsion subplots in.
:type fig: `matplotlib.figure.Figure`
:param tors1: List of torsions from ligand 1.
:type tors1: list[fep_edge_data.FEPTorsions]
:param tors2: List of torsions from ligand 2.
:type tors2: list[fep_edge_data.FEPTorsions]
:param tors_from: Starting torsion index.
:type tors_from: int
:param tors_to: Ending torsion index
:type tors_to: int
:param ipage: Page index of the torsions, only needed for interactively
highlighting torsion atoms.
:type ipage: int or None
:param structure1_item: Structure item 1, only needed for interactively
highlighting torsion atoms.
:type structure1_item: `structure2d.structure_item` or None
:param structure2_item: Structure item 2, only needed for interactively
highlighting torsion atoms.
:type structure2_item: `structure2d.structure_item` or None
:param for_print: Whether the figure is used for print.
:type for_print: bool
:returns: A list of the energy line plot axes
:rtype: list[matplotlib.axes.Axes]
"""
fig.clear()
gs = gridspec.GridSpec(tors_to - tors_from,
2,
top=0.90,
bottom=0.1,
left=0.05,
right=0.95,
wspace=0.1,
hspace=0.3,
width_ratios=[3., 3.])
tors1_page = []
tors2_page = []
torsion_colors = []
line_plots = []
for irow in range(tors_from, tors_to):
tor1, tor2 = None, None
if irow < len(tors1):
tor1 = tors1[irow]
tors1_page.append(tor1)
lig1_ax = fig.add_subplot(gs[(irow % TORSIONS_PER_PAGE) * 2])
if irow < len(tors2):
tor2 = tors2[irow]
tors2_page.append(tor2)
if not tor1:
lig1_ax = None
lig2_ax = fig.add_subplot(gs[(irow % TORSIONS_PER_PAGE) * 2 +
1],
sharey=lig1_ax)
color = rh.colors[irow % len(rh.colors)]
torsion_colors.append(QtGui.QColor(color))
hist_tick_pos = None
if irow == tors_from:
hist_tick_pos = 'top'
if tor1:
if irow in [tors_to - 1, len(tors1) - 1]:
hist_tick_pos = 'bottom'
tor1_line_plot = tor1.plot(lig1_ax,
color=color,
for_print=for_print,
pot_tick_pos='left',
hist_tick_pos=hist_tick_pos)
line_plots.append(tor1_line_plot)
if tor2:
if irow in [tors_to - 1, len(tors2) - 1]:
hist_tick_pos = 'bottom'
tor2_line_plot = tor2.plot(lig2_ax,
color=color,
for_print=for_print,
pot_tick_pos='right',
hist_tick_pos=hist_tick_pos)
line_plots.append(tor2_line_plot)
if ipage is not None:
self.gs_dict[ipage] = gs
self.fig_dict[ipage] = fig
self.lig1_item_dict[ipage] = structure1_item
self.lig2_item_dict[ipage] = structure2_item
fig.canvas.mpl_connect('motion_notify_event', self.on_move)
self.tors1_dict[ipage] = tors1_page
self.tors2_dict[ipage] = tors2_page
self.torsion_colors_dict[ipage] = torsion_colors
return line_plots
def _get_torsion_plot_img(self, fig, nrows):
fig.set_size_inches(7.25, 0.65 * nrows)
temp_plot = self._get_temp_image_fn()
fig.savefig(temp_plot, bbox_inches='tight', dpi=300)
size = fig.get_size_inches()
(sx, sy) = rh.aspectScale(size[0], size[1], 7.25, 7.0)
img = rh.platypus.Image(temp_plot, sx * rh.inch, sy * rh.inch)
return img
def _ligand_torsion_report(self):
"""
Generates report which monitors Rotatable bonds
"""
Ele = self.Elements
fed = self.fed
lig1_tors, lig2_tors = fed.ligand_torsions
torsion_figs = []
lig2d_annotated = []
tors_per_page = []
npages = 1
if lig1_tors.tors_total > 0 or lig2_tors.tors_total > 0:
maxtors = max(lig1_tors.tors_total, lig2_tors.tors_total)
npages += maxtors // TORSIONS_PER_PAGE
if maxtors % (npages * TORSIONS_PER_PAGE) == TORSIONS_PER_PAGE:
npages -= 1
for ipage in range(npages):
tors_from = ipage * TORSIONS_PER_PAGE
tors_to = min(maxtors, (ipage + 1) * TORSIONS_PER_PAGE)
tors_fig = self.get_torsions_plot(lig1_tors.all_tors,
lig2_tors.all_tors,
tors_from,
tors_to,
for_print=True)
torsion_figs.append(tors_fig)
lig1_2d_annotated, lig2_2d_annotated = \
self.get_2d_tors_annotated_lig_pair(lig1_tors.all_tors,
lig2_tors.all_tors,
tors_from,
tors_to)
lig1_2d_img = self._get_annotated_2d_lig_img(lig1_2d_annotated)
lig2_2d_img = self._get_annotated_2d_lig_img(lig2_2d_annotated)
lig2d_annotated.append((lig1_2d_img, lig2_2d_img))
tors_per_page.append(tors_to - tors_from)
for ipage in range(npages):
rh.new_page(Ele)
rh.report_add_logo(Ele)
if ipage == 0:
rh.header(Ele, "Ligand Conformation Analysis")
else:
rh.header(Ele,
"Ligand Conformation Analysis (Cont. %i)" % ipage)
# add table with annotated 2d ligands
ligs_2d_image = lig2d_annotated[ipage]
table2 = [[ligs_2d_image[0], ligs_2d_image[1]]]
table2.append(['Ligand 1', 'Ligand 2'])
width2 = [3.5, 3.5]
# yapf: disable
style2 = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 1), (-1, -1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table2, style2, width2)
try:
Ele.append(
self._get_torsion_plot_img(torsion_figs[ipage],
tors_per_page[ipage]))
except:
rh.pargph(Ele, 'Torsion profiles cannot be generated')
rh.pargph(Ele, self._get_ligand_torsion_text())
def _get_ligand_torsion_text(self):
if self.fed.perturbation_type == FEP_TYPES.LIGAND_SELECTIVITY:
return FEPEdgeReportMaker._get_ligand_torsion_text()
text = """\
Rotatable bonds (Rb) in both ligands are enumerated and color-coded. For
each Rb, a representative dihedral angle is monitored throughout the
complex and solvent simulation legs, their distributions are then
plotted. Hollow bars show <b>solvent</b> and filled bars show
<b>complex</b> leg distributions. Input starting conformation is marked
as a <font color="gray">gray</font> vertical line.
Potential energy around each Rb overlays the plot with the dark-blue
curve and corresponding labels on the Y-axis. Local strain energies are
shown below the plot. The units are in <i>kcal/mol</i>.
"""
return text
def _ligand_report(self):
"""
Generate ligand report
"""
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Ligand Analysis for End-Point λ-Replicas")
if fed.ligand1_charge > 0:
str_lig1_charge = '+' + str(fed.ligand1_charge)
else:
str_lig1_charge = str(fed.ligand1_charge)
if fed.ligand2_charge > 0:
str_lig2_charge = '+' + str(fed.ligand2_charge)
else:
str_lig2_charge = str(fed.ligand2_charge)
str_lig1_mass = "%.1f au" % fed.ligand1_atomic_mass
str_lig2_mass = "%.1f au" % fed.ligand2_atomic_mass
table = [[
' ', 'PDB Name', 'No. Atoms', 'No. Heavy', 'No. Hot Atoms',
'Rot. Bonds', 'Atomic Mass', 'Charge'
]]
table.append([
'Ligand 1:', fed.ligand1_pdb_name, fed.ligand1_total_atoms,
fed.ligand1_total_heavy, fed.ligand1_total_hot,
fed.ligand1_rot_bonds, str_lig1_mass, str_lig1_charge
])
table.append([
'Ligand 2:', fed.ligand2_pdb_name, fed.ligand2_total_atoms,
fed.ligand2_total_heavy, fed.ligand2_total_hot,
fed.ligand2_rot_bonds, str_lig2_mass, str_lig2_charge
])
width = [0.8] + 7 * [0.93]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('ALIGN', (0, 1), (0, nrows - 1), 'RIGHT'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
self._ligand_rmsd_wrt_protein()
self._ligand_properties()
def _ligand_properties(self):
Ele = self.Elements
fed = self.fed
rh.add_spacer(Ele)
rh.pargph(Ele, '<u>Ligand Misc. Properties</u>')
rh.add_spacer(Ele)
def wrap(mean_stdev):
return rh.get_pargph("%.1f ± %.2f" % mean_stdev,
fontsize=9,
hAlign='center')
def plain_wrap(text):
return rh.get_pargph(r"%s" % text, fontsize=9, hAlign='center')
table = [['', '', 'Ligand 1', '', 'Ligand2', ''],
[
'', 'Units', 'In solution', 'In complex', 'In solution',
'In complex'
]]
table.append([
'RMSD',
plain_wrap('Å'),
wrap(fed.ligand1_sol_sid_rmsd()),
wrap(fed.ligand1_cpx_sid_rmsd()),
wrap(fed.ligand2_sol_sid_rmsd()),
wrap(fed.ligand2_cpx_sid_rmsd())
])
table.append([
'Radius of gyration',
plain_wrap('Å'),
wrap(fed.ligand1_sol_sid_rgyr()),
wrap(fed.ligand1_cpx_sid_rgyr()),
wrap(fed.ligand2_sol_sid_rgyr()),
wrap(fed.ligand2_cpx_sid_rgyr())
])
table.append([
'Molecular SA',
plain_wrap('Å<sup>2</sup>'),
wrap(fed.ligand1_sol_sid_molsa()),
wrap(fed.ligand1_cpx_sid_molsa()),
wrap(fed.ligand2_sol_sid_molsa()),
wrap(fed.ligand2_cpx_sid_molsa())
])
table.append([
'Solvent-accessible SA',
plain_wrap('Å<sup>2</sup>'),
wrap(fed.ligand1_sol_sid_sasa()),
wrap(fed.ligand1_cpx_sid_sasa()),
wrap(fed.ligand2_sol_sid_sasa()),
wrap(fed.ligand2_cpx_sid_sasa())
])
table.append([
'Polar SA',
plain_wrap('Å<sup>2</sup>'),
wrap(fed.ligand1_sol_sid_psa()),
wrap(fed.ligand1_cpx_sid_psa()),
wrap(fed.ligand2_sol_sid_psa()),
wrap(fed.ligand2_cpx_sid_psa())
])
table.append([
'Intramolecular HB',
plain_wrap('#'),
wrap(fed.ligand1_sol_sid_lighb()),
wrap(fed.ligand1_cpx_sid_lighb()),
wrap(fed.ligand2_sol_sid_lighb()),
wrap(fed.ligand2_cpx_sid_lighb())
])
table.append([
'Number of waters',
plain_wrap('#'),
plain_wrap("N/A"),
wrap(fed.ligand1_cpx_sid_waters()),
plain_wrap("N/A"),
wrap(fed.ligand2_cpx_sid_waters())
])
table.append([
'Ligand strain',
plain_wrap('kcal/mol'),
wrap(fed.ligand1_sol_sid_rb_strain()),
wrap(fed.ligand1_cpx_sid_rb_strain()),
wrap(fed.ligand2_sol_sid_rb_strain()),
wrap(fed.ligand2_cpx_sid_rb_strain())
])
ncol = len(table[0])
nrows = len(table)
width = [1.6, 0.8] + 4 * [1.0]
# yapf: disable
style = [('SPAN', (2, 0), (3, 0)),
('SPAN', (4, 0), (5, 0)),
('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (1, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (ncol - 1, 1), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
def _ligand_rmsd_wrt_protein(self, name1='Ligand 1', name2='Ligand 2'):
Ele = self.Elements
fed = self.fed
rh.pargph(Ele, '<u>Ligand RMSD in complex</u>')
rmsd1 = fed.receptor_sid_rmsd_ligand_lambda0
rmsd2 = fed.receptor_sid_rmsd_ligand_lambda1
rmsd_plot = self._get_ligand_wrt_prot_rmsd_plot(
rmsd1,
rmsd2,
fed.cpx_sid_snapshot_times_ps,
fed.cpx_sim_time,
name1=name1,
name2=name2)
rmsd_text = self._get_ligand_wrt_prot_rmsd_text()
table = [[rmsd_plot, rh.get_pargph(rmsd_text, fontsize=10)]]
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('VALIGN', (0, 0), (-1, -1), 'TOP')]
# yapf: enable
width = [4.3, 3.]
rh.add_vtable(Ele, table, style, width)
def _get_ligand_wrt_prot_rmsd_text(self):
text = """\
The Root Mean Square Deviation (RMSD) of a ligand is measured with
respect to the input ligand pose by first aligning the complex on the
protein. If the values observed are significantly larger than the RMSD
of the protein (See the <i>Protein Analysis Report</i>), then it is
likely that the ligand has diffused away from its initial binding site.
Remember that since the FEP calculation is in the <i>'REST'</i> mode, the
fluctuations may be larger than observed for typical <i>MD</i> jobs.
"""
return text
def _get_ligand_wrt_prot_rmsd_plot(self, rmsd_l0, rmsd_l1, ts, ts_max,
name1, name2):
names = [name1, name2]
ylabel = r"Ligand RMSD ($\AA$)"
y_height = 2.0
img = self._get_pair_rmsd_plot(rmsd_l0,
rmsd_l1,
ts,
ts_max,
names,
ylabel,
y_height,
for_print=True)
return img
def _protein_report(self):
"""
Generates the protein report (RSMD/RSMF for both end-point lambdas).
"""
Ele = self.Elements
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Protein Analysis for End-Point λ-Replicas")
self._protein_details()
rh.add_spacer(Ele)
try:
Ele.append(self._get_sequence_viewer_image())
except:
rh.pargph(Ele, '<i>Protein sequence could not be generated.</i>')
self._protein_rmsd()
self._protein_rmsf()
def _protein_rmsf(self):
Ele = self.Elements
rh.pargph(Ele, '<u>Protein RMSF</u>')
rh.add_spacer(Ele)
show_common_contacts = True
show_uniq_lig1 = True
show_uniq_lig2 = True
show_sse = True
show_b_factor = self._b_factor_valid(self.fed.receptor_b_factor)
contacts = (show_common_contacts, show_uniq_lig1, show_uniq_lig2)
rmsf_plot = self.get_protein_rmsf(for_print=True,
show_sse=show_sse,
show_b_factor=show_b_factor,
show_interacting_residues=contacts)
Ele.append(rmsf_plot)
rh.pargph(Ele,
self._get_protein_rmsf_text(show_sse=show_sse,
show_b=show_b_factor,
lig_res=contacts),
fontsize=10)
def _get_protein_rmsf_text(self,
show_sse=False,
show_b=False,
lig_res=None):
"""
Returns a string to be used in the PDF document
"""
text = """ \
The Root Mean Square Fluctuation (RMSF) is useful for characterizing
local changes along the protein chain. Protein backbone RMSF is shown
in this plot. Typically you will observe that the tails (<i>N</i>- and
<i>C</i>-terminal) fluctuate more than any other part of the protein.
Secondary structure elements like alpha helices and beta strands are
usually more rigid than the unstructured part of the protein, and thus
fluctuate less than the loop regions.
"""
if show_b:
text += """\
Experimental <font color='#F781F3'>B factors</font> extracted from
PDB is overlaid alongside RMSF values. RMSF and B factors
fluctuations should correlate, but not necessarily follow each
other.
"""
if show_sse:
text += """ \
<br/><font color='#F06040'>Alpha-helical</font> and
<font color='#40C0E0'>beta-strand</font> regions are highlighted in
red and blue backgrounds, respectively. These regions are defined
by helices or strands that persist over 70% of the entire
simulation. One should expect secondary structure elements to be in
the low-fluctuating RMSF region.
"""
if lig_res:
text += """ \
Vertical lines indicate the residues that are involved in
interacting with the ligand.
"""
return text
def _add_sse_span(self, ax, data_limits, color):
"""
Setup overlay for SSE elements
:param ax: matplotlib plot
:type ax: matplotlib.ax
:param data_limits: list of ranges for secondary structure elements
along the residue index
:type data_limits: `(from,to)`
:param color: which color to use
:type color: matplotlib color, either string or hex
"""
for (i, j) in data_limits:
ax.axvspan(i, j, facecolor=color, alpha=0.1)
def _add_b_factor(self, ax2, b_data, color, y_range_ratio, fsize=8):
"""
Add B factor to the plot, on the Y2 axis. scale the Y2 axis to be
compatible with Y1 axis.
:param ax2: matplotlib plot
:type ax2: matplotlib.ax (twinx/Y2 axis)
:param b_data: b factor values for each residue
:type b_data: `list` of `float`
:param color: which color to use
:type color: matplotlib color, either string or hex
:param y_range_ratio: a ratio by which to scale the Y2-axis
:type y_range_ratio: float
:param fsize: font size to use in plots
:type fsize: int
:return: a title of the plot for legends
:rtype: str
"""
p, = ax2.plot(range(len(b_data)), b_data, color=color, label='B factor')
ax2.set_ylabel('B factor',
rotation=270,
color=color,
size=fsize,
alpha=0.9)
b_aver = numpy.median(b_data)
b_ymax = b_aver * (y_range_ratio - 1)
ax2.set_ylim(ymin=b_ymax / -10, ymax=b_ymax)
# do not show zero tick
ax2.yaxis.set_major_locator(MaxNLocator(nbins=7, prune='lower'))
# color each tick magenta
[i.set_color(color) for i in ax2.get_yticklabels()]
return p
def _add_interacting_residues(self, ax, show_res, vals, color, ymin=0):
"""
Draw a vertical bar for reisdues that interact with the ligand
:param ax: matplotlib plot
:type ax: matplotlib.ax
:param show_res: residues index to add the vbar to
:type show_res: int
:param vals: the hight of the vbar
:type vals: float
:param color: which color to use
:type color: matplotlib color, either string or hex
:param ymin: the starting y-position of the bar
:type ymin: float
"""
ax.vlines(show_res, ymin, vals, colors=color, linestyle='solid')
[docs] def gen_protein_rmsf_image(self,
rmsf_l0,
rmsf_l1,
sse_l0,
sse_l1,
b_factor,
contact_data,
show_res_type=None):
fig = figure.Figure(figsize=(7.25, 2.25))
FigureCanvas(fig)
ax = fig.add_subplot(111)
self._gen_protein_rmsf_ax(ax,
rmsf_l0,
rmsf_l1,
sse_l0,
sse_l1,
b_factor,
contact_data,
show_res_type,
fsize=8)
rh.change_plot_colors(ax)
temp_plot = self._get_temp_image_fn()
fig.savefig(temp_plot, bbox_inches='tight', dpi=300)
size = fig.get_size_inches()
(sx, sy) = rh.aspectScale(size[0], size[1], 7.25, 2.25)
img = rh.platypus.Image(temp_plot, sx * rh.inch, sy * rh.inch)
return img
def _get_rmsf_contacts(self, show_interacting_residues):
"""
Return the residue contact type strings, which are either 'common',
'uniq_lig1' and 'uniq_lig2' as well as a dictionary using those strings
as keys that maps them to the list of residue indices for the contact
type.
:param show_interacting_residues: which types of interactions (common,
unique to ligand1, and unique to ligand2) that need to be obtained
:type show_interacting_residues: three-tuple of bool
:return: the lig res type list and the contact data dictionary
:rtype: tuple(list(str), dict(str, list))
"""
res_contact_data = {}
if True not in show_interacting_residues:
return [], {}
protein_residues1, protein_residues2 = self.fed.get_protein_residues_sequence_tags(
)
lig1_res = self.fed.receptor_residues_interaction_ligand1
lig2_res = self.fed.receptor_residues_interaction_ligand2
lig1_idx = set([protein_residues1.index(r) for r in lig1_res])
lig2_idx = set([protein_residues2.index(r) for r in lig2_res])
common = lig1_idx.intersection(lig2_idx)
res_contact_data['common'] = list(common)
res_contact_data['uniq_lig1'] = list(
lig1_idx.symmetric_difference(common))
res_contact_data['uniq_lig2'] = list(
lig2_idx.symmetric_difference(common))
types = ['common', 'uniq_lig1', 'uniq_lig2']
lig_res_type = []
for show, inter_type in zip(show_interacting_residues, types):
if show:
lig_res_type.append(inter_type)
return lig_res_type, res_contact_data
[docs] def get_protein_rmsf(self,
for_print=True,
fig=None,
show_sse=False,
show_b_factor=False,
show_interacting_residues=(False, False, False)):
"""
This function returns either an img or a matplotlib object of protein
RMSF plot
:param for_print: return an image for PDF if True otherwise return
matplotlib object
:type for_print: bool
:param fig: matplotlib Figure object
:type fig: `matplotlib.figure.Figure`
:param show_sse: overlay secondary structure information on RMSF plot
:type show_sse: bool
:param show_b_factor: overlay b factor on the y2 axis
:type show_b_factor: bool
:param show_interacting_residues: a tuple with residue indices that
correspond to different contacts for ligand1, ligand2, and common
:type show_interacting_residues: (common, uniq_lig1, uniq_lig2)
:return: returns either an image or an matplotlib object, depends on the
for_print variable
:rtype: image | matplotlib ax obj
"""
lambda0_rmsf = self.fed.receptor_sid_rmsf_backbone_lambda0
lambda1_rmsf = self.fed.receptor_sid_rmsf_backbone_lambda1
lambda0_sse = {self.HELIX: None, self.STRAND: None}
lambda1_sse = {self.HELIX: None, self.STRAND: None}
if show_sse:
helix, strand = self.fed.sse_limits_lambda0
lambda0_sse = {self.HELIX: helix, self.STRAND: strand}
helix, strand = self.fed.sse_limits_lambda1
lambda1_sse = {self.HELIX: helix, self.STRAND: strand}
b_factor = None
if show_b_factor:
b_factor = self.fed.receptor_b_factor
lig_res_type, res_contact_data = self._get_rmsf_contacts(
show_interacting_residues)
if for_print:
return self.gen_protein_rmsf_image(lambda0_rmsf,
lambda1_rmsf,
lambda0_sse,
lambda1_sse,
b_factor,
res_contact_data,
show_res_type=lig_res_type)
else:
if fig is None:
fig = figure.Figure()
ax = fig.add_subplot(111)
else:
fig.canvas.ax.clear()
ax = fig.canvas.ax
# Clear twin axis (B-factor) and set its visibility.
if hasattr(ax, 'ax2'):
ax.ax2.clear()
ax.ax2.get_yaxis().set_visible(show_b_factor)
self._gen_protein_rmsf_ax(ax,
lambda0_rmsf,
lambda1_rmsf,
lambda0_sse,
lambda1_sse,
b_factor,
res_contact_data,
show_res_type=lig_res_type,
fsize=12)
return ax
def _b_factor_valid(self, b_factor: List[float]):
"""
If b factor values are not available don't plot them
"""
if not b_factor:
return False
if max(b_factor) - min(b_factor) < 1:
return False
return True
def _gen_protein_rmsf_ax(self,
ax,
rmsf_l0,
rmsf_l1,
sse_l0,
sse_l1,
b_factor,
contact_data,
show_res_type=None,
fsize=8):
"""
Given an matplotlib ax object, setup the plot
:param rmsf_l0: list of rmsf valuse for protein in lambda0 sim
:type rmsf_l0: `float`
:param rmsf_l1: list of rmsf valuse for protein in lambda1 sim
:type rmsf_l1: `float`
:param sse_l0: A dictionary with 'helix' and 'strand' values which
contain limits of all the SSE in lambda0 protein
:type sse_l0: dict{'helix_limits', 'strand_limits'}
:param sse_l1: A dictionary with 'helix' and 'strand' values which
contain limits of all the SSE in lambda1 protein
:type sse_l1: dict{'helix_limits', 'strand_limits'}
:param b_factor: b factor values for each residue
:type b_factor: `list` of `float`
:param contact_data: a dictionary with residue indices that correspond
to different contacts for ligand1, ligand2, and common
:type contact_data: dict{'common', 'uniq_lig1', 'uniq_lig2'}
:param show_res_type: to show contact residues
:type show_res_type: bool
"""
rmsf_val = [rmsf_l0, rmsf_l1]
names = [r'$\lambda$ = 0', r'$\lambda$ = 1']
colors = ['black', 'brown']
lines = []
tmp_list = rmsf_l0 + rmsf_l1
maxval = max(tmp_list)
aver = numpy.mean(tmp_list)
ylim = maxval
if maxval > aver * 4:
ylim = aver * 4
res_index = list(range(len(rmsf_l0)))
for data, name, color in zip(rmsf_val, names, colors):
p, = ax.plot(res_index, data, color=color, label=name, alpha=0.9)
lines.append(p)
if sse_l0[self.HELIX] is not None and sse_l1[self.HELIX] is not None:
self._add_sse_span(ax, sse_l0[self.HELIX] + sse_l1[self.HELIX],
self.HELIX_COLOR)
if sse_l0[self.STRAND] is not None and sse_l1[self.STRAND] is not None:
self._add_sse_span(ax, sse_l0[self.STRAND] + sse_l1[self.STRAND],
self.STRAND_COLOR)
if b_factor is not None:
if hasattr(ax, 'ax2'):
ax2 = ax.ax2
else:
ax2 = ax.twinx()
ax.ax2 = ax2
y_range_ratio = numpy.max(rmsf_val) / numpy.mean(rmsf_val)
p = self._add_b_factor(ax2, b_factor, self.B_FACTOR, y_range_ratio,
fsize)
lines.append(p)
ax2.tick_params(axis='y', colors=self.B_FACTOR, labelsize=fsize)
if show_res_type and contact_data:
colors = {
'common': 'red',
'uniq_lig1': 'black',
'uniq_lig2': 'brown'
}
for kw in show_res_type:
vals = []
for i in contact_data[kw]:
vals.append(max(rmsf_l0[i - 1], rmsf_l1[i - 1]))
self._add_interacting_residues(ax,
contact_data[kw],
vals,
colors[kw],
ymin=ylim / -10)
ax.legend(lines,
names,
loc=4,
borderaxespad=0.6,
ncol=2,
borderpad=0.2,
handlelength=1,
handletextpad=0.2,
columnspacing=0.4,
fancybox=True,
shadow=False,
prop={'size': fsize})
ax.set_ylim(ymin=ylim / -10, ymax=ylim)
ax.set_xlim(xmin=0, xmax=len(res_index))
leg = ax.get_legend()
[item.set_linewidth(4) for item in leg.legendHandles]
ax.set_xlabel("Residue Index", size=fsize)
ax.set_ylabel(r"Protein RMSF ($\AA$)", size=fsize)
ax.yaxis.set_major_locator(MaxNLocator(nbins=7, prune='lower'))
def _protein_rmsd(self):
Ele = self.Elements
fed = self.fed
rh.add_spacer(Ele)
rh.pargph(Ele, '<u>Protein RMSD</u>')
lambda0 = self.fed.receptor_sid_rmsd_backbone_lambda0
lambda1 = self.fed.receptor_sid_rmsd_backbone_lambda1
rmsd_plot = self._get_protein_lambda_rmsd_plot(
lambda0, lambda1, fed.cpx_sid_snapshot_times_ps, fed.cpx_sim_time)
rmsd_text = self._get_protein_rmsd_text()
table = [[rmsd_plot, rh.get_pargph(rmsd_text, fontsize=10)]]
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('VALIGN', (0, 0), (-1, -1), 'TOP')]
# yapf: enable
width = [4.3, 3.]
rh.add_vtable(Ele, table, style, width)
def _get_protein_rmsd_text(self):
text = """ \
The Root Mean Square Deviation (RMSD) of a protein is measured here
over the backbone atoms. Monitoring protein RMSD for two end-points in
the FEP simulation may give insights into its structural stability.
Values of the order of 1-3 Å are perfectly acceptable for
medium-size, globular proteins. Changes much larger than that, however,
indicate that the protein is undergoing a larger than expected
conformational changes for an equilibrated system, during the
simulation.
"""
return text
def _get_protein_lambda_rmsd_plot(self, rmsd_l0, rmsd_l1, ts, ts_max):
names = [r'$\lambda$ = 0', r'$\lambda$ = 1']
ylabel = r"Protein RMSD ($\AA$)"
img = self._get_pair_rmsd_plot(rmsd_l0,
rmsd_l1,
ts,
ts_max,
names,
ylabel,
for_print=True)
return img
def _get_pair_rmsd_plot(self,
val0,
val1,
ts,
ts_max,
names,
ylabel=r'RMSD ($\AA$)',
y_inch=1.8,
for_print=False):
"""
:param val0: list of RSMD values in angstroms for lambda=0
:type val0: `float`
:param val1: list of RSMD values in angstroms for lambda=1
:type val1: `float`
:param ts: list of simulation times for each snapshot
:type ts: `float`
:param ts_max: simulation time
:type ts_max: float
:param names: List of names that that are in the legend
:type names: `str`
:param ylabel: Label that goes on the y-axis
:type ylabel: str
:param y_inch: height of the image returned
:type y_inch: float
:param for_print: Bool if the figure will be used for print
:type for_print: bool
"""
if val0 is not None and val1 is not None and len(val0) != len(val1):
raise InconsistentDataException(
"Inconsistent number of RMSD values")
fig = figure.Figure(figsize=(4.2, y_inch))
FigureCanvas(fig)
self._create_pair_rmsd_fig(val0,
val1,
ts,
ts_max,
names,
fig,
ylabel,
for_print=for_print)
temp_plot = self._get_temp_image_fn()
fig.savefig(temp_plot, bbox_inches='tight', dpi=300)
size = fig.get_size_inches()
(sx, sy) = rh.aspectScale(size[0], size[1], 4.2, y_inch)
img = rh.platypus.Image(temp_plot, sx * rh.inch, sy * rh.inch)
return img
def _create_pair_rmsd_fig(self,
val0,
val1,
ts,
ts_max,
names,
fig,
ylabel=r'RMSD ($\AA$)',
for_print=False):
fontsize = 12
if for_print:
fontsize = 8
fig.clear()
ax = fig.add_subplot(111)
vals = [val0, val1]
colors = ['black', 'brown']
lines = []
for data, name, color in zip(vals, names, colors):
if data is None:
continue
p, = ax.plot(ts, data, color=color, label=name, alpha=0.9)
lines.append(p)
ax.legend(lines,
names,
loc=4,
borderaxespad=0.6,
ncol=2,
borderpad=0.2,
handlelength=1,
handletextpad=0.2,
columnspacing=0.4,
fancybox=True,
shadow=False,
prop={'size': fontsize})
ax.set_ylim(ymin=0)
ax.set_xlim(xmax=max(ts_max, ts[-1]))
leg = ax.get_legend()
[item.set_linewidth(4) for item in leg.legendHandles]
ax.set_xlabel("Time (nsec)", size=fontsize)
ax.set_ylabel(ylabel, size=fontsize)
ax.yaxis.set_major_locator(MaxNLocator(nbins=7, prune='lower'))
if for_print:
rh.change_plot_colors(ax)
def _get_sequence_viewer_image(self):
"""
Generate protein sequence image, given a protein structure
Returns a platypus image
"""
seq_viewer = self.SeqViewer()
seq_viewer.incorporateStructure(self.fed.receptor_st)
seq_viewer.wrapped = True
seq_viewer.has_ruler = False
seq_viewer.display_identity = False
seq_viewer.display_boundaries = True
seq_viewer.chain_name_only = True
seq_viewer.hide_empty_lines = True
seq_viewer.crop_image = True
seq_viewer.wrapped_width = 90
seq_viewer.setColorMode(self.COLOR_POSITION)
seq_viewer.addAnnotation(self.ANNOTATION_SSBOND, remove=True)
seq_viewer.addAnnotation(self.ANNOTATION_RESNUM)
temp_seq_file = self._get_temp_image_fn()
size = seq_viewer.saveImage(temp_seq_file)
(sx, sy) = rh.aspectScale(size[0], size[1], 7.0, 5.0)
protein_img = rh.platypus.Image(temp_seq_file, sx * rh.inch,
sy * rh.inch)
return protein_img
def _free_energy_convergence_profiles(self, table=True):
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Free Energy Convergence")
rh.pargph(Ele, '<u>%s Leg</u>' % fed.leg1_name)
rh.add_spacer(Ele)
sol_convergence = self._gen_free_energy_convergence_plot(
fed.sol_start_time_ns, fed.sol_end_time_ns, fed.sol_delta_g_forward,
fed.sol_delta_g_forward_err, fed.sol_delta_g_reverse,
fed.sol_delta_g_reverse_err, fed.sol_delta_g_sliding,
fed.sol_delta_g_sliding_err)
Ele.append(sol_convergence)
if table:
self._gen_free_energy_convergence_table(
fed.sol_delta_g_forward_df_per_replica,
fed.sol_delta_g_forward_dg,
fed.sol_delta_g_forward_bootstrap_std,
fed.sol_delta_g_forward_analytical_std) # yapf: disable
rh.add_spacer(Ele)
rh.add_spacer(Ele)
rh.pargph(Ele, '<u>%s Leg</u>' % fed.leg2_name)
rh.add_spacer(Ele)
cpx_convergence = self._gen_free_energy_convergence_plot(
fed.cpx_start_time_ns, fed.cpx_end_time_ns, fed.cpx_delta_g_forward,
fed.cpx_delta_g_forward_err, fed.cpx_delta_g_reverse,
fed.cpx_delta_g_reverse_err, fed.cpx_delta_g_sliding,
fed.cpx_delta_g_sliding_err)
Ele.append(cpx_convergence)
if table:
self._gen_free_energy_convergence_table(
fed.cpx_delta_g_forward_df_per_replica,
fed.cpx_delta_g_forward_dg,
fed.cpx_delta_g_forward_bootstrap_std,
fed.cpx_delta_g_forward_analytical_std) # yapf: disable
rh.add_spacer(Ele)
rh.pargph(Ele, self._get_free_energy_convergence_text(table=table))
def _gen_free_energy_convergence_table(self, df_per_replica, dg,
bootstrap_std, analytical_std):
nreplicas = len(df_per_replica) + 1
if nreplicas <= 12:
fontsize, col_width, head_width = 11, 0.50, 1.2
elif nreplicas <= 18:
fontsize, col_width, head_width = 11, 0.42, 1.2
else:
fontsize, col_width, head_width = 6.5, 0.3, 0.7
replica_name = [
rh.get_pargph('λ pair:', color='#888888', fontsize=fontsize)
]
replica_dF = ['Forward dF:']
replica_dE1 = ['Bootstrap STD:']
replica_dE2 = ['Analytical STD:']
for i, val in enumerate(df_per_replica):
name = "%i,%i" % (i, i + 1)
replica_name.append(name)
replica_dF.append('%.2f' % float(val[0]))
replica_dE1.append('%.3f' % float(val[1]))
replica_dE2.append('%.3f' % float(val[2]))
# add total dG
replica_name.append('Total')
replica_dF.append('%.2f' % dg)
replica_dE1.append('%.3f' % bootstrap_std)
replica_dE2.append('%.3f' % analytical_std)
color = rh.rlcolors.green if dg < 0 else rh.rlcolors.orange
table = [replica_name, replica_dF, replica_dE1, replica_dE2]
width = [head_width] + (nreplicas * [col_width])
nfields = len(width)
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('FONTSIZE', (0, 0), (-1, -1), fontsize),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('ALIGN', (0, 0), (0, nrows - 1), 'LEFT'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray),
('BOX', (-1, 0), (-1, -1), 0.4, color)]
# yapf: enable
rh.add_vtable(self.Elements, table, style, width)
def _ligand_perturbation_details(self):
Ele = self.Elements
fed = self.fed
if fed.ligand1_charge > 0:
str_lig1_charge = '+' + str(fed.ligand1_charge)
else:
str_lig1_charge = str(fed.ligand1_charge)
if fed.ligand2_charge > 0:
str_lig2_charge = '+' + str(fed.ligand2_charge)
else:
str_lig2_charge = str(fed.ligand2_charge)
str_lig1_mass = "%.1f au" % fed.ligand1_atomic_mass
str_lig2_mass = "%.1f au" % fed.ligand2_atomic_mass
rh.pargph(Ele, '<u>Ligand Information</u>')
table = [[
' ', 'Name', 'HexID', 'No. Atoms / Heavy', 'Atomic Mass', 'Charge',
'Mol. Formula'
]]
table.append([
'Ligand 1:', fed.ligand1_name, fed.ligand1_hash,
"%i / %i" % (fed.ligand1_total_atoms, fed.ligand1_total_heavy),
str_lig1_mass, str_lig1_charge, fed.ligand1_mol_formula
])
table.append([
'Ligand 2:', fed.ligand2_name, fed.ligand2_hash,
"%i / %i" % (fed.ligand2_total_atoms, fed.ligand2_total_heavy),
str_lig2_mass, str_lig2_charge, fed.ligand2_mol_formula
])
width = [0.8, 1.6, 0.6, 1.4, 0.8, 0.6, 1.0]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('ALIGN', (0, 1), (0, nrows - 1), 'RIGHT'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.add_spacer(Ele)
lig1_2d_image, lig2_2d_image = self._generate_2d_lig_images()
table2 = [['Ligand1', 'Ligand2']]
table2.append([lig1_2d_image, lig2_2d_image])
width2 = [3.5, 3.5]
nfields2 = len(table2[0])
# yapf: disable
style2 = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields2 - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table2, style2, width2)
def _generate_2d_lig_images(self):
"""
Generate 2d images of the ligands. Try to align them first, or just
regular 2d images, if nothing works, return a string.
"""
try:
lig1_2d_image, lig2_2d_image = self._generate_aligned_2d_lig_pair()
except:
try:
lig1_2d_image = self.get_2d_ligand_image(self.fed.ligand1_st)
lig2_2d_image = self.get_2d_ligand_image(self.fed.ligand2_st)
except:
lig1_2d_image = "Couldn't generate 2d plot"
lig2_2d_image = "Couldn't generate 2d plot"
return lig1_2d_image, lig2_2d_image
def _generate_aligned_2d_lig_pair(self):
"""
Generate *aligned* 2d images of the ligands.
"""
fed = self.fed
fn = [
self._get_temp_image_fn('lig1_'),
self._get_temp_image_fn('lig2_')
]
rh.generate_aligned_2d_ligand_pair(fn, fed.ligand1_st, fed.ligand2_st,
fed._atom_mapping)
img1 = rh.get_image(fn[0], 3.25, 3.25)
img2 = rh.get_image(fn[1], 3.25, 3.25)
return img1, img2
def _protein_details(self):
Ele = self.Elements
fed = self.fed
rh.pargph(Ele, '<u>Protein Information</u>')
rh.add_spacer(Ele)
charge = fed.receptor_charge
if charge > 0:
charge = '+' + str(charge)
table = [[
'Name', 'Tot. Residues', 'Prot. Chain(s)', 'Res. in Chain(s)',
'No. Atoms', 'No. Heavy Atoms', 'Charge'
]]
# receptor title may be too long to fit in the pdf column
receptor_title = fed.receptor_title if len(fed.receptor_title) < 20 \
else f'{fed.receptor_title[:20]}...'
table.append([
receptor_title, fed.receptor_total_residues,
fed.receptor_chain_names, fed.receptor_total_residues_in_chains,
fed.receptor_total_atom, fed.receptor_total_heavy, charge
])
nfields = len(table[0])
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, 1.1)
def _salt_details(self):
Ele = self.Elements
fed = self.fed
if fed.sol_salt_info is None or fed.cpx_salt_info is None:
return
rh.add_spacer(Ele)
rh.pargph(Ele, '<u>Salt Information</u>')
rh.add_spacer(Ele)
table = [[
'Ion Type', 'Concentration in Complex', 'Concentration in Solvent'
]]
sol = {ion[0]: ion[1] for ion in fed.sol_salt_info}
cpx = {ion[0]: ion[1] for ion in fed.cpx_salt_info}
ions = set(list(sol.keys()) + list(cpx.keys()))
for ion in ions:
table.append((ion, sol.get(ion, '0.0 mM'), cpx.get(ion, '0.0 mM')))
nrows = len(table)
nfields = len(table[0])
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, [1., 2., 2.])
def _report_fep_simulation_details(self):
Ele = self.Elements
fed = self.fed
rh.pargph(Ele, '<u>System Details</u>')
rh.add_spacer(Ele)
rh.pargph(Ele,
"<font color='#888888'>Job name</font>: %s" % fed.jobname)
rh.pargph(
Ele, "<font color='#888888'>Job type</font>: %s" %
fed.perturbation_type)
if self.perturbation_type in [
FEP_TYPES.SMALL_MOLECULE, FEP_TYPES.COVALENT_LIGAND,
FEP_TYPES.METALLOPROTEIN, FEP_TYPES.ABSOLUTE_BINDING
]:
rh.pargph(
Ele,
"<font color='#888888'>Perturbation</font>: (%s ↔ %s)" %
(fed.ligand1_name, fed.ligand2_name))
else:
rh.pargph(
Ele, "<font color='#888888'>Perturbation</font>: %s" %
(fed.prm_name))
rh.add_spacer(Ele)
sol_temp_str = '%0.1f' % fed.sol_temperature
cpx_temp_str = '%0.1f' % fed.cpx_temperature
sol_atoms_water = str(fed.sol_total_atoms) + ' / ' + str(
fed.sol_total_waters)
cpx_atoms_water = str(fed.cpx_total_atoms) + ' / ' + str(
fed.cpx_total_waters)
sol_delta_g = rh.get_pargph("%.2f ± %.2f" % fed.sol_delta_g,
fontsize=10)
cpx_delta_g = rh.get_pargph("%.2f ± %.2f" % fed.cpx_delta_g,
fontsize=10)
delta_g_text = rh.get_pargph("ΔG [kcal/mol]",
fontsize=10,
color='gray')
table = [[
'', 'Ensemble', 'Force Field', 'Temp. [K]', 'Sim Time [ns]',
'No. Atoms / Waters', delta_g_text
]]
table.append([
'%s:' % fed.leg1_name, fed.sol_ensemble, fed.sol_ff, sol_temp_str,
fed.sol_sim_time, sol_atoms_water, sol_delta_g
])
table.append([
'%s:' % fed.leg2_name, fed.cpx_ensemble, fed.cpx_ff, cpx_temp_str,
fed.cpx_sim_time, cpx_atoms_water, cpx_delta_g
])
if fed.ddg_corrections_map:
correction = rh.get_pargph(
"%.2f ± %.2f" %
(fed.ddg_corrections_sum.val, fed.ddg_corrections_sum.unc),
fontsize=10)
table.append(['Corr. Term:', '', '', '', '', '', correction])
width = [1.4, 0.8, 0.7, 0.7, 0.9, 1.3, 1.1]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (1, 0), (-1, -1), 'CENTER'),
('ALIGN', (0, 1), (0, nrows - 1), 'RIGHT'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.header(
Ele, "Difference in binding free energy (ΔΔG) "
"is: %.2f kcal/mol" % fed.delta_delta_g)
[docs] def error(self, msg):
print(str(msg))
[docs] def str(self, str_in):
"""
This is to remove the "'s in ARK returned strings enclosed in
double quotes
"""
return str(str_in).strip("\"")
[docs]class PRMEdgeReportMaker(FEPEdgeReportMaker):
"""
This class generates a PDF report for FEP+ for protein residue mutation.
"""
[docs] def __init__(self,
fep_edge_data,
basename=None,
perturbation_type=FEP_TYPES.SMALL_MOLECULE):
"""
This class generates a PDF report for an FEP/REST (+) type of job. Both
Legs of the simulations are processed in the `FEPEdgeData` and is used
by this class to compile a result.
:type fep_edge_data: `FEPEdgeData`
:param fep_edge_data: Object containing all the data for this report
:type basename: string
:param basename: the basename of the file of the PDF report
:type perturbation_type: str
:param perturbation_type: Type of PRM job it was
"""
FEPEdgeReportMaker.__init__(self, fep_edge_data, basename,
perturbation_type)
[docs] def report(self):
self._init_report()
self._report_fep_simulation_details()
if self.perturbation_type == FEP_TYPES.LIGAND_SELECTIVITY:
self._ligand_details()
self._fragment_perturbation_details()
self._protein_details()
self._salt_details()
self._free_energy_convergence_profiles()
self._replica_exchange_density()
self._protein_report()
if self.perturbation_type not in [
FEP_TYPES.PROTEIN_STABILITY, FEP_TYPES.PROTEIN_SELECTIVITY
]:
self._ligand_report()
self._ligand_torsion_report()
try:
self._protein_ligand_report()
except:
pass
self.doc.build(self.Elements, canvasmaker=rh.NumberedCanvas)
self._cleanup_temp_files()
print(self.basename, 'report is written')
def _ligand_details(self):
Ele = self.Elements
fed = self.fed
str_lig_mass = "%.1f au" % fed.ligand_mass
str_lig_charge = "%i" % fed.ligand_charge
if fed.ligand_charge > 0:
str_lig_charge = "+%i" % fed.ligand_st.formal_charge
rh.pargph(Ele, '<u>Ligand Information</u>')
table = [[
' ', 'Name', 'No. Atoms / Heavy', 'Atomic Mass', 'Charge',
'Mol. Formula'
]]
table.append([
'Ligand:', fed.ligand_name,
"%i / %i" % (fed.ligand_st.atom_total, fed.ligand_nheavy),
str_lig_mass, str_lig_charge, fed.ligand_st_mol_formula
])
width = [1.6, 1.0, 1.4, 1.0, 1.0]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('ALIGN', (0, 1), (0, nrows - 1), 'RIGHT'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.add_spacer(Ele)
lig_2d_image = self.get_2d_ligand_image(self.fed.ligand_st)
table_lig = [['Ligand']]
table_lig.append([lig_2d_image])
width_lig = [3.5]
nfields_lig = len(table_lig[0])
# yapf: disable
style2 = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields_lig - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table_lig, style2, width_lig)
def _fragment_perturbation_details(self):
Ele = self.Elements
fed = self.fed
if fed.ligand1_charge > 0:
str_lig1_charge = '+' + str(fed.ligand1_charge)
else:
str_lig1_charge = str(fed.ligand1_charge)
if fed.ligand2_charge > 0:
str_lig2_charge = '+' + str(fed.ligand2_charge)
else:
str_lig2_charge = str(fed.ligand2_charge)
str_lig1_mass = "%.1f au" % fed.ligand1_atomic_mass
str_lig2_mass = "%.1f au" % fed.ligand2_atomic_mass
rh.pargph(Ele, '<u>Peptide Fragment Information</u>')
table = [[
' ', 'Name', 'HexID', 'No. Atoms / Heavy', 'Atomic Mass', 'Charge',
'Mol. Formula'
]]
table.append([
'Frag. 1:', fed.wt_name, fed.ligand1_hash,
"%i / %i" % (fed.ligand1_total_atoms, fed.ligand1_total_heavy),
str_lig1_mass, str_lig1_charge, fed.ligand1_mol_formula
])
table.append([
'Frag. 2:', fed.mut_name, fed.ligand2_hash,
"%i / %i" % (fed.ligand2_total_atoms, fed.ligand2_total_heavy),
str_lig2_mass, str_lig2_charge, fed.ligand2_mol_formula
])
width = [0.8, 2.0, 0.6, 1.4, 0.8, 0.6, 1.0]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('ALIGN', (0, 1), (0, nrows - 1), 'RIGHT'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.add_spacer(Ele)
lig1_2d_image, lig2_2d_image = self._generate_2d_lig_images()
table2 = [['Fragment 1', 'Fragment 2']]
table2.append([lig1_2d_image, lig2_2d_image])
width2 = [3.5, 3.5]
nfields2 = len(table2[0])
# yapf: disable
style2 = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields2 - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table2, style2, width2)
def _protein_ligand_report(self):
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Fragment Interactions for End-Point "
"λ-Replicas")
residues = fed.cpx_sid_protein_residues
pli0 = fed.cpx_sid_pl_results[0]
pli1 = fed.cpx_sid_pl_results[1]
nframes = fed.cpx_sid_number_of_frames
st1, st2 = fed.wt_frag_st, fed.mut_frag_st
name1, name2 = fed.wt_name, fed.mut_name
if self.perturbation_type not in [
FEP_TYPES.PROTEIN_STABILITY, FEP_TYPES.PROTEIN_SELECTIVITY
]:
st1, st2 = fed.ligand_st, fed.ligand_st
name1, name2 = 'Ligand', 'Ligand'
bar_chart = self._get_plot_pl_bar_chart(residues,
pli0,
pli1,
nframes,
name_1=name1,
name_2=name2)
Ele.append(bar_chart)
rh.pargph(Ele, self._get_pl_text(), fontsize=10)
lp_stats = fed.cpx_sid_lp_results
cutoff = 0.20
lig_noh_atom_dict, lig_atom_dict = fed._get_ligand_atom_dict(st1)
atom_mapping = fed._frag_atom_mapping
lid0_img = self._get_lid_img(lp_stats[0],
lig_noh_atom_dict,
lig_atom_dict,
st1,
atom_mapping[0],
template_st=None,
template_core_idx=None,
cutoff=cutoff)
lig_noh_atom_dict, lig_atom_dict = fed._get_ligand_atom_dict(st2)
lid1_img = self._get_lid_img(lp_stats[1],
lig_noh_atom_dict,
lig_atom_dict,
st2,
atom_mapping[1],
template_st=st1,
template_core_idx=atom_mapping[0],
cutoff=cutoff)
table = [[name1, name2], [lid0_img, lid1_img]]
width = [4.0, 4.0]
nfields = len(table[0])
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
def _get_pl_text(self):
text = FEPEdgeReportMaker._get_pl_text(self)
if self.perturbation_type in [
FEP_TYPES.PROTEIN_STABILITY, FEP_TYPES.PROTEIN_SELECTIVITY
]:
return text.replace('ligand', 'residue')
return text.replace('ligand', 'fragment')
def _get_free_energy_convergence_text(self, table=True):
text = FEPEdgeReportMaker._get_free_energy_convergence_text(self, table)
if self.perturbation_type in [
FEP_TYPES.PROTEIN_STABILITY, FEP_TYPES.PROTEIN_SELECTIVITY
]:
return text.replace('ligand', 'residue')
return text
def _get_ligand_torsion_text(self):
return """\
Rotatable bonds (Rb) of a ligand in each end point lambda windows are
enumerated and color-coded. For each Rb, a representative dihedral angle
is monitored throughout the simulation for both complex and solvent
legs. The distributions of these conformations is then plotted for each
Rb. In addition, potential energy around each Rb overlays the histograms
with the dark-blue curve and corresponding labels on the <i>Y</i>-axis.
The energy units are reported in <i>kcal/mol</i>.
"""
def _getTorsionStructures(self):
if self.perturbation_type == FEP_TYPES.LIGAND_SELECTIVITY:
return self.fed.ligand_st, self.fed.ligand_st
return self.fed.wt_frag_st, self.fed.mut_frag_st
def _generate_2d_lig_images(self):
"""
Generate 2d images of the ligands. Try to align them first, or just
regular 2d images, if nothing works, return a string.
"""
try:
lig1_2d_image, lig2_2d_image = self._generate_aligned_2d_lig_pair()
except:
try:
st1, st2 = self.fed.wt_frag_st, self.fed.mut_frag_st
lig1_2d_image = self.get_2d_ligand_image(st1)
lig2_2d_image = self.get_2d_ligand_image(st2)
except:
lig1_2d_image = "Couldn't generate 2d plot"
lig2_2d_image = "Couldn't generate 2d plot"
return lig1_2d_image, lig2_2d_image
def _generate_aligned_2d_lig_pair(self):
"""
Generate *aligned* 2d images of the ligands.
"""
fed = self.fed
fn = [
self._get_temp_image_fn('lig1_'),
self._get_temp_image_fn('lig2_')
]
st1, st2 = self.fed.wt_frag_st, self.fed.mut_frag_st
rh.generate_aligned_2d_ligand_pair(fn, st1, st2, fed._atom_mapping)
img1 = rh.get_image(fn[0], 3.25, 3.25)
img2 = rh.get_image(fn[1], 3.25, 3.25)
return img1, img2
def _ligand_properties(self):
Ele = self.Elements
fed = self.fed
rh.add_spacer(Ele)
rh.pargph(Ele, '<u>Ligand Misc. Properties</u>')
rh.add_spacer(Ele)
def wrap(mean_stdev):
return rh.get_pargph("%.1f ± %.2f" % mean_stdev,
fontsize=9,
hAlign='center')
def plain_wrap(text):
return rh.get_pargph(r"%s" % text, fontsize=9, hAlign='center')
table = [['', 'Units', 'In WT', 'In Mutant']]
table.append([
'RMSD',
plain_wrap('Å'),
wrap(fed.ligand1_cpx_sid_rmsd()),
wrap(fed.ligand2_cpx_sid_rmsd())
])
table.append([
'Radius of gyration',
plain_wrap('Å'),
wrap(fed.ligand1_cpx_sid_rgyr()),
wrap(fed.ligand2_cpx_sid_rgyr())
])
table.append([
'Molecular SA',
plain_wrap('Å<sup>2</sup>'),
wrap(fed.ligand1_cpx_sid_molsa()),
wrap(fed.ligand2_cpx_sid_molsa())
])
table.append([
'Solvent-accessible SA',
plain_wrap('Å<sup>2</sup>'),
wrap(fed.ligand1_cpx_sid_sasa()),
wrap(fed.ligand2_cpx_sid_sasa())
])
table.append([
'Polar SA',
plain_wrap('Å<sup>2</sup>'),
wrap(fed.ligand1_cpx_sid_psa()),
wrap(fed.ligand2_cpx_sid_psa())
])
table.append([
'Intramolecular HB',
plain_wrap('#'),
wrap(fed.ligand1_cpx_sid_lighb()),
wrap(fed.ligand2_cpx_sid_lighb())
])
table.append([
'Number of waters',
plain_wrap('#'),
wrap(fed.ligand1_cpx_sid_waters()),
wrap(fed.ligand2_cpx_sid_waters())
])
table.append([
'Ligand strain',
plain_wrap('kcal/mol'),
wrap(fed.ligand1_cpx_sid_rb_strain()),
wrap(fed.ligand2_cpx_sid_rb_strain())
])
ncol = len(table[0])
nrows = len(table)
width = [1.6, 0.8] + 4 * [1.0]
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (1, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (ncol - 1, 1), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
def _protein_details(self):
Ele = self.Elements
fed = self.fed
rh.pargph(Ele, '<u>Protein Information</u>')
rh.add_spacer(Ele)
charge = fed.receptor_charge
if charge > 0:
charge = '+' + str(charge)
table = [[
'Name', 'Tot. Residues', 'Prot. Chain(s)', 'Res. in Chain(s)',
'No. Atoms', 'No. Heavy Atoms', 'Charge'
]]
table.append([
fed.receptor_title, fed.receptor_total_residues,
fed.receptor_chain_names, fed.receptor_total_residues_in_chains,
fed.receptor_total_atom, fed.receptor_total_heavy, charge
])
nfields = len(table[0])
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, 1.1)
def _ligand_report(self):
"""
Generate ligand report
"""
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Ligand Analysis for End-Point λ-Replicas")
rh.add_spacer(Ele)
str_lig_charge = str(fed.ligand_charge)
if fed.ligand_charge > 0:
str_lig_charge = '+' + str_lig_charge
str_lig_mass = "%.1f au" % fed.ligand_mass
table = [[
' ', 'PDB Name', 'No. Atoms', 'No. Heavy', 'Rot. Bonds',
'Atomic Mass', 'Charge'
]]
table.append([
'Ligand:', fed.ligand_name, fed.ligand_st.atom_total,
fed.ligand_nheavy, fed.ligand_rb_total, str_lig_mass, str_lig_charge
])
width = [0.8] + 7 * [0.93]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('ALIGN', (0, 1), (0, nrows - 1), 'RIGHT'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.add_spacer(Ele)
self._ligand_rmsd_wrt_protein(name1='In WT', name2='In Mutant')
self._ligand_properties()
[docs]class ABFEPReportMaker(FEPEdgeReportMaker):
[docs] def report(self):
self._init_report()
self._report_fep_simulation_details()
self._ligand_perturbation_details()
self._protein_details()
self._salt_details()
self._free_energy_convergence_profiles(table=False)
self._replica_exchange_density()
self._protein_report()
self._ligand_report()
self._ligand_torsion_report()
self._protein_ligand_report()
self.doc.build(self.Elements, canvasmaker=rh.NumberedCanvas)
self._cleanup_temp_files()
print(self.basename, 'report is written')
def _report_fep_simulation_details(self):
"""
Report primary information about the simulation
"""
Ele = self.Elements
fed = self.fed
rh.pargph(Ele, '<u>System Details</u>')
rh.add_spacer(Ele)
rh.pargph(Ele,
"<font color='#888888'>Job name</font>: %s" % fed.jobname)
rh.pargph(
Ele, "<font color='#888888'>Job type</font>: %s" %
fed.perturbation_type)
rh.pargph(
Ele, "<font color='#888888'>Perturbation</font>: (%s ↔ %s)" %
(fed.ligand1_name, 'None'))
rh.add_spacer(Ele)
sol_temp_str = '%0.1f' % fed.sol_temperature
cpx_temp_str = '%0.1f' % fed.cpx_temperature
sol_atoms_water = f'{fed.sol_total_atoms} / {fed.sol_total_waters}'
cpx_atoms_water = f'{fed.cpx_total_atoms} / {fed.cpx_total_waters}'
sol_delta_g = rh.get_pargph("%.2f ± %.2f" % fed.sol_delta_g,
fontsize=10)
cpx_delta_g = rh.get_pargph("%.2f ± %.2f" % fed.cpx_delta_g,
fontsize=10)
delta_g_text = rh.get_pargph("ΔG [kcal/mol]",
fontsize=10,
color='gray')
table = [[
'', 'Ensemble', 'Force Field', 'Temp. [K]', 'Sim Time [ns]',
'No. Atoms / Waters', delta_g_text
]]
table.append([
'%s:' % fed.leg1_name, fed.sol_ensemble, fed.sol_ff, sol_temp_str,
fed.sol_sim_time, sol_atoms_water, sol_delta_g
])
table.append([
'%s:' % fed.leg2_name, fed.cpx_ensemble, fed.cpx_ff, cpx_temp_str,
fed.cpx_sim_time, cpx_atoms_water, cpx_delta_g
])
if fed.ddg_corrections_map:
correction = rh.get_pargph(
"%.2f ± %.2f" %
(fed.ddg_corrections_sum.val, fed.ddg_corrections_sum.unc),
fontsize=10)
table.append(['Corr. Term:', '', '', '', '', '', correction])
width = [1.4, 0.8, 0.7, 0.7, 0.9, 1.3, 1.1]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (1, 0), (-1, -1), 'CENTER'),
('ALIGN', (0, 1), (0, nrows - 1), 'RIGHT'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.header(
Ele, "Absolute binding free energy (ΔΔG) "
"is: %.2f kcal/mol" % fed.delta_delta_g)
def _ligand_perturbation_details(self):
Ele = self.Elements
fed = self.fed
rh.pargph(Ele, '<u>Ligand Information</u>')
table = [[
'Name', 'HexID', 'No. Atoms / Heavy', 'Atomic Mass', 'Charge',
'Mol. Formula'
]]
table.append([
fed.ligand1_name, fed.ligand1_hash,
f'{fed.ligand1_total_atoms} / {fed.ligand1_total_heavy}',
f'{fed.ligand1_atomic_mass:.1f} au',
f'{"+" if fed.ligand1_charge > 0 else ""}{fed.ligand1_charge}',
fed.ligand1_mol_formula
])
width = [2.2, 0.8, 1.6, 1.2, 0.8, 1.5]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.add_spacer(Ele)
try:
lig1_2d_image = self.get_2d_ligand_image(self.fed.ligand1_st)
except:
lig1_2d_image = "Couldn't generate 2d plot"
table2 = [['Ligand'], [lig1_2d_image]]
width2 = [4.0]
nfields2 = len(table2[0])
# yapf: disable
style2 = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields2 - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table2, style2, width2)
def _protein_rmsf(self):
Ele = self.Elements
rh.pargph(Ele, '<u>Protein RMSF</u>')
rh.add_spacer(Ele)
show_sse = True
show_b_factor = self._b_factor_valid(self.fed.receptor_b_factor)
rmsf_plot = self.get_protein_rmsf(for_print=True,
show_sse=show_sse,
show_b_factor=show_b_factor,
show_interacting_residues=(False,
False,
False))
Ele.append(rmsf_plot)
rh.pargph(Ele,
self._get_protein_rmsf_text(show_sse=show_sse,
show_b=show_b_factor,
lig_res=False),
fontsize=10)
def _ligand_report(self):
"""
Generate ligand report
"""
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Ligand Analysis for End-Point λ-Replicas")
table = [[
'Title ', 'PDB Name', 'No. Atoms', 'No. Heavy', 'No. Hot Atoms',
'Rot. Bonds', 'Atomic Mass', 'Charge'
]]
table.append([
fed.ligand1_name, fed.ligand1_pdb_name, fed.ligand1_total_atoms,
fed.ligand1_total_heavy, fed.ligand1_total_hot,
fed.ligand1_rot_bonds, f'{fed.ligand1_atomic_mass:.1f} au',
f'{"+" if fed.ligand1_charge > 0 else ""}{fed.ligand1_charge}'
])
width = [1.5] + 7 * [0.85]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.add_spacer(Ele)
self._ligand_rmsd_wrt_protein()
self._ligand_properties()
def _ligand_rmsd_wrt_protein(self, name1='Ligand'):
Ele = self.Elements
fed = self.fed
rh.pargph(Ele, '<u>Ligand RMSD in complex</u>')
rmsd_plot = self._get_ligand_wrt_prot_rmsd_plot(
fed.receptor_sid_rmsd_ligand_lambda0,
None,
fed.cpx_sid_snapshot_times_ps,
fed.cpx_sim_time,
name1=name1,
name2=None)
rmsd_text = self._get_ligand_wrt_prot_rmsd_text()
table = [[rmsd_plot, rh.get_pargph(rmsd_text, fontsize=10)]]
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('VALIGN', (0, 0), (-1, -1), 'TOP')]
# yapf: enable
width = [4.3, 3.]
rh.add_vtable(Ele, table, style, width)
def _ligand_properties(self):
Ele = self.Elements
fed = self.fed
rh.add_spacer(Ele)
rh.pargph(Ele, '<u>Ligand Misc. Properties</u>')
rh.add_spacer(Ele)
def wrap(mean_stdev):
return rh.get_pargph("%.1f ± %.2f" % mean_stdev,
fontsize=9,
hAlign='center')
def plain_wrap(text):
return rh.get_pargph(r"%s" % text, fontsize=9, hAlign='center')
table = [['', 'Units', 'Solvent Leg', 'Complex Leg']]
table.append([
'RMSD',
plain_wrap('Å'),
wrap(fed.ligand1_sol_sid_rmsd()),
wrap(fed.ligand1_cpx_sid_rmsd()),
])
table.append([
'Radius of gyration',
plain_wrap('Å'),
wrap(fed.ligand1_sol_sid_rgyr()),
wrap(fed.ligand1_cpx_sid_rgyr()),
])
table.append([
'Molecular SA',
plain_wrap('Å<sup>2</sup>'),
wrap(fed.ligand1_sol_sid_molsa()),
wrap(fed.ligand1_cpx_sid_molsa()),
])
table.append([
'Solvent-accessible SA',
plain_wrap('Å<sup>2</sup>'),
wrap(fed.ligand1_sol_sid_sasa()),
wrap(fed.ligand1_cpx_sid_sasa()),
])
table.append([
'Polar SA',
plain_wrap('Å<sup>2</sup>'),
wrap(fed.ligand1_sol_sid_psa()),
wrap(fed.ligand1_cpx_sid_psa()),
])
table.append([
'Intramolecular HB',
plain_wrap('#'),
wrap(fed.ligand1_sol_sid_lighb()),
wrap(fed.ligand1_cpx_sid_lighb()),
])
table.append([
'Number of waters',
plain_wrap('#'),
plain_wrap("N/A"),
wrap(fed.ligand1_cpx_sid_waters()),
])
table.append([
'Ligand strain',
plain_wrap('kcal/mol'),
wrap(fed.ligand1_sol_sid_rb_strain()),
wrap(fed.ligand1_cpx_sid_rb_strain()),
])
ncol = len(table[0])
nrows = len(table)
width = [1.6, 0.8] + 2 * [1.0]
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (1, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (ncol - 1, 1), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
def _ligand_torsion_report(self):
"""
Generates report which monitors Rotatable bonds
"""
TORSIONS_PER_PAGE = 10
Ele = self.Elements
fed = self.fed
lig1_tors, lig2_tors = fed.ligand_torsions
torsion_figs = []
lig2d_annotated = []
tors_per_page = []
npages = 1
if lig1_tors.tors_total > 0:
ntors = lig1_tors.tors_total
npages += lig1_tors.tors_total // TORSIONS_PER_PAGE
if ntors % (npages * TORSIONS_PER_PAGE) == TORSIONS_PER_PAGE:
npages -= 1
for ipage in range(npages):
tors_from = ipage * TORSIONS_PER_PAGE
tors_to = min(ntors, (ipage + 1) * TORSIONS_PER_PAGE)
tors_fig = self.get_torsions_plot(lig1_tors.all_tors, [],
tors_from,
tors_to,
for_print=True)
torsion_figs.append(tors_fig)
lig1_2d_annotated, lig2_2d_annotated = \
self.get_2d_tors_annotated_lig_pair(lig1_tors.all_tors,
lig2_tors.all_tors,
tors_from,
tors_to)
lig1_2d_img = self._get_annotated_2d_lig_img(lig1_2d_annotated)
lig2d_annotated.append(lig1_2d_img)
tors_per_page.append(tors_to - tors_from)
for ipage in range(npages):
rh.new_page(Ele)
rh.report_add_logo(Ele)
if ipage == 0:
rh.header(Ele, "Ligand Conformation Analysis")
else:
rh.header(Ele,
"Ligand Conformation Analysis (Cont. %i)" % ipage)
# add table with annotated 2d ligands
table2 = [[lig2d_annotated[ipage]], [fed.ligand1_name]]
width2 = [4.0]
# yapf: disable
style2 = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 1), (-1, -1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table2, style2, width2)
try:
Ele.append(
self._get_torsion_plot_img(torsion_figs[ipage],
tors_per_page[ipage]))
except:
rh.pargph(Ele, 'Torsion profiles cannot be generated')
rh.pargph(Ele, self._get_ligand_torsion_text())
[docs] def create_torsions_plot(self,
fig,
tors1,
tors2,
tors_from,
tors_to,
ipage=None,
structure1_item=None,
structure2_item=None,
for_print=False):
"""
Creates a plot for torsions using the given matplot figure.
:param fig: Figure to draw the torsion subplots in.
:type fig: `matplotlib.figure.Figure`
:param tors1: List of torsions from ligand 1.
:type tors1: list[fep_edge_data.FEPTorsions]
:param tors2: List of torsions from ligand 2.
:type tors2: list[fep_edge_data.FEPTorsions]
:param tors_from: Starting torsion index.
:type tors_from: int
:param tors_to: Ending torsion index
:type tors_to: int
:param ipage: Page index of the torsions, only needed for interactively
highlighting torsion atoms.
:type ipage: int or None
:param structure1_item: Structure item 1, only needed for interactively
highlighting torsion atoms.
:type structure1_item: `structure2d.structure_item` or None
:param structure2_item: Structure item 2, only needed for interactively
highlighting torsion atoms.
:type structure2_item: `structure2d.structure_item` or None
:param for_print: Whether the figure is used for print.
:type for_print: bool
:returns: A list of the energy line plot axes
:rtype: list[matplotlib.axes.Axes]
"""
fig.clear()
gs = gridspec.GridSpec(math.ceil((tors_to - tors_from) / 2),
2,
top=0.90,
bottom=0.1,
left=0.05,
right=0.95,
wspace=0.1,
hspace=0.3,
width_ratios=[3., 3.])
tors_page = []
torsion_colors = []
line_plots = []
for itor in range(tors_from, tors_to):
tor = tors1[itor]
tors_page.append(tor)
lig_ax = fig.add_subplot(gs[(itor % TORSIONS_PER_PAGE)])
color = rh.colors[itor % len(rh.colors)]
torsion_colors.append(QtGui.QColor(color))
hist_tick_pos = None
if itor in [tors_from, tors_from + 1]:
hist_tick_pos = 'top'
# show tick labels location for the last set of dihedrals
if (itor == tors_to - 1) or \
(not ((tors_to - tors_from) % 2) and itor == tors_to - 2):
hist_tick_pos = 'bottom'
tor_line_plot = tor.plot(lig_ax,
color=color,
for_print=for_print,
pot_tick_pos='right' if itor %
2 else 'left',
hist_tick_pos=hist_tick_pos)
line_plots.append(tor_line_plot)
if ipage is not None:
self.gs_dict[ipage] = gs
self.fig_dict[ipage] = fig
self.lig1_item_dict[ipage] = structure1_item
self.lig2_item_dict[ipage] = structure2_item
fig.canvas.mpl_connect('motion_notify_event', self.on_move)
self.tors1_dict[ipage] = tors_page
self.tors2_dict[ipage] = tors_page
self.torsion_colors_dict[ipage] = torsion_colors
return line_plots
[docs] def on_move(self, event):
# override the super class because we now have two columns of torsion
# plots, both of which are associated with one ligand.
page = self.torsions_page_index
lig1_item = self.lig1_item_dict[page]
if event.xdata is None or event.ydata is None:
# outside the subplots
self.remove_last_atom_highlight(lig1_item)
return
# get the grid positions (bottom, top, left and right locations)
# relative to the entire figure, but converted into pixels.
fig = self.fig_dict[page]
fig_width, fig_height = fig.get_size_inches() * fig.dpi
bots, tops, lefts, rights = self.gs_dict[page].get_grid_positions(fig)
bots, tops = bots * fig_height, tops * fig_height
lefts, rights = lefts * fig_width, rights * fig_width
# The number of rows and columns per page
nrow, ncol = len(bots), len(lefts)
for row_idx in range(nrow):
if bots[row_idx] <= event.y <= tops[row_idx]:
# cursor is over a row that has torsion plots.
# now determine the index in the torsion dictionary that the
# cursor is in
for col_idx in range(ncol):
if lefts[col_idx] <= event.x <= rights[col_idx]:
break
idx = 2 * row_idx + col_idx
color = self.torsion_colors_dict[page][idx]
tors1 = self.tors1_dict[page]
self._highlight_matching_torsion_atoms(idx, tors1, lig1_item,
color)
return
def _get_rmsf_contacts(self, show_interacting_residues):
# override the super class because we only need to honor the 'uniq_lig1'
# (index of 1) request which will show all contacts with the one ligand.
# There is no concept of 'common' nor 'uniq_lig2' interactions since
# there is no second ligand of interest.
if not show_interacting_residues[1]:
return [], {}
protein_residues = self.fed.receptor_residue_sequence_tags
contact_residues = self.fed.receptor_residues_interaction_ligand1
contacts = [protein_residues.index(r) for r in contact_residues]
return ['uniq_lig1'], {'uniq_lig1': contacts}
def _protein_ligand_report(self):
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Protein-Ligand Interactions")
residues = fed.cpx_sid_protein_residues
pli0, pli1 = fed.cpx_sid_pl_results[0], []
nframes = fed.cpx_sid_number_of_frames
bar_chart = self._get_plot_pl_bar_chart(residues, pli0, pli1, nframes)
Ele.append(bar_chart)
rh.pargph(Ele, self._get_pl_text(), fontsize=10)
rh.add_spacer(Ele)
lp_stats = fed.cpx_sid_lp_results
cutoff = 0.20
lig_noh_atom_dict, lig_atom_dict = fed.get_ligand1_atom_dict()
lid0_img = self._get_lid_img(lp_stats[0],
lig_noh_atom_dict,
lig_atom_dict,
fed.ligand1_st, [],
template_st=None,
template_core_idx=None,
cutoff=cutoff)
table = [['Ligand'], [lid0_img]]
width = [4.0]
nfields = len(table[0])
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
def _create_pl_bar_fig(self,
label,
data_dict0,
data_dict1,
nframes,
pl_interactions=ALL_PL_INTERACTIONS,
name_1="Ligand 1",
name_2="Ligand 2",
fig=None,
for_print=False):
nresidues = len(label)
rrange = numpy.arange(nresidues)
sum_list0 = numpy.zeros(nresidues)
if fig is None:
fig = figure.Figure(figsize=(7, 3))
else:
fig.clear()
ax = fig.add_subplot(111)
for type_ in pl_interactions:
itype_vals0 = numpy.array([data_dict0[lab][type_] for lab in label
]) / float(nframes) * 100
ax.bar(rrange,
itype_vals0,
bottom=sum_list0,
label=type_,
color=PL_INTERACTIONS[type_][0])
sum_list0 += itype_vals0
fontsize = 8 if for_print else 11
#ax1.yaxis.tick_right()
ax.xaxis.set_ticks(rrange)
# Check if residue labels have the same chain name
new_label = [lab[2:].replace('_', ' ') for lab in label] \
if len(set(lab[0] for lab in label)) == 1 \
else [lab.replace('_', ' ') for lab in label]
ax.xaxis.set_ticklabels(new_label, rotation=70, fontsize=fontsize)
ax.set_xlim(-0.5, nresidues - 0.5)
ax.set_ylim(0, math.ceil(max(sum_list0) / 100) * 100)
ax.set_ylabel('% interact. with Ligand', fontsize=fontsize)
if for_print:
for tick in ax.get_xticklabels():
tick.set_fontsize(7.5)
rh.change_plot_colors(ax)
return fig
def _get_pl_text(self):
text = """\
Above bar chart illustrates the type of interactions the protein
residues make with the ligand. Multiple types of specific interactions
are monitored throughout the simulation and provide a way to examine and
compare how the protein interacts with the ligand.
The specific interactions types monitored and displayed are:
<b><font color='#7FC97F'> hydrogen bond</font></b>,
<b><font color='#BEAED4'>hydrophobic</font></b>,
<b><font color='#F0027F'>ionic</font></b> and
<b><font color='#386CB0'>water bridges</font></b>. The stacked bar charts
are normalized over the course of the trajectory. More information
about the geometry and the different interaction subtype categories
can be found in the Schrödinger's <i>Desmond User Manual</i>,
under 'Simulation Interactions Diagram' (SID) section. <i><b>Note:</b></i>
The values may exceed 100% as the residue can make multiple contacts of
the same type with the ligand.
"""
return text
[docs]class SolubilityReportMaker(FEPEdgeReportMakerBase):
'''
This class generates plots, images and PDF reports for Solubility FEP jobs.
'''
VERBOSE = False
TORSIONS_PER_PAGE = 10
[docs] def __init__(self,
fep_edge_data,
basename=None,
perturbation_type=FEP_TYPES.SOLUBILITY):
"""
This class generates a PDF report for an FEP/REST (+) type of job. Both
Legs of the simulations are processed in the `FEPEdgeData` and is used
by this class to compile a result.
:type fep_edge_data: `FEPEdgeData`
:param fep_edge_data: Object containing all the data for this report
:type basename: string
:param basename: the basename of the file of the PDF report
:type perturbation_type: str
:param perturbation_type: FEP Type
"""
super(SolubilityReportMaker, self).__init__(fep_edge_data, basename,
perturbation_type)
self._load_modules()
self.fig_dict = {}
self.gs_dict = {}
self.mol_item_dict = {}
self.hydration_tors_dict = {}
self.sublimation_tors_dict = {}
self.torsion_colors_dict = {}
self.torsions_page_index = 0
[docs] def report(self):
self._init_report()
self._report_fep_simulation_details()
self._molecule_perturbation_details()
self._free_energy_convergence_profiles()
self._replica_exchange_density()
self._molecule_report()
self._molecule_environment_report()
self._mol_torsion_report()
self.doc.build(self.Elements, canvasmaker=rh.NumberedCanvas)
self._cleanup_temp_files()
print(self.basename, 'report is written')
def _init_report(self):
pdf_filename = self.basename + '.pdf'
self.doc = rh.platypus.SimpleDocTemplate(pdf_filename)
self.doc.title = u'FEP+ Solubility Report by Schrodinger Inc.'
self.doc.author = "Dmitry Lupyan, PhD"
self.doc.leftMargin = 30.
self.doc.rightMargin = 20.
self.doc.topMargin = 10.
self.doc.bottomMargin = 10.
rh.report_add_logo(self.Elements)
rh.header(self.Elements,
"FEP+ Solubility Simulation Details and Results")
def _mol_torsion_report(self):
"""
Generates report which monitors Rotatable bonds
"""
Ele = self.Elements
fed = self.fed
mol_tors = fed.mol_torsions
torsion_figs = []
mol2d_annotated = []
tors_per_page = []
npages = 1
if mol_tors[0].tors_total > 0:
ntors = mol_tors[0].tors_total
npages += ntors // self.TORSIONS_PER_PAGE
if ntors % (npages *
self.TORSIONS_PER_PAGE) == self.TORSIONS_PER_PAGE:
npages -= 1
tors_figs = []
for ipage in range(npages):
tors_from = ipage * self.TORSIONS_PER_PAGE
tors_to = min(ntors, (ipage + 1) * self.TORSIONS_PER_PAGE)
tors_figs = self.get_torsions_plot(mol_tors,
tors_from,
tors_to,
for_print=True)
torsion_figs.append(tors_figs)
tors_per_page.append(tors_to - tors_from)
mol_2d_annotated, _ = \
self.get_2d_tors_annotated_lig_pair(mol_tors[0].all_tors,
mol_tors[0].all_tors,
tors_from,
tors_to)
mol_2d_img = self._get_annotated_2d_lig_img(mol_2d_annotated)
mol2d_annotated.append(mol_2d_img)
for ipage in range(npages):
rh.new_page(Ele)
rh.report_add_logo(Ele)
if ipage == 0:
rh.header(Ele, "Ligand Conformation Analysis")
else:
rh.header(Ele, f"Ligand Conformation Analysis (Cont. {ipage})")
# add table with annotated 2d ligands
table2 = [[mol2d_annotated[ipage]], [fed.mol_name]]
width2 = [4.0]
# yapf: disable
style2 = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 1), (-1, -1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table2, style2, width2)
try:
Ele.append(
self._get_torsion_plot_img(torsion_figs[ipage],
tors_per_page[ipage]))
except:
rh.pargph(Ele, 'Torsion profiles cannot be generated')
rh.pargph(Ele, self._get_mol_torsion_text())
def _get_torsion_plot_img(self, fig, nrows):
fig.set_size_inches(7.25, 0.65 * nrows)
temp_plot = self._get_temp_image_fn()
fig.savefig(temp_plot, bbox_inches='tight', dpi=300)
size = fig.get_size_inches()
(sx, sy) = rh.aspectScale(size[0], size[1], 7.25, 7.0)
img = rh.platypus.Image(temp_plot, sx * rh.inch, sy * rh.inch)
return img
[docs] def get_torsions_plot(self, all_tors, tors_from, tors_to, for_print=False):
fig = figure.Figure()
FigureCanvas(fig)
# This is for generating PDF report, no interaction needed.
self.create_torsions_plot(fig,
all_tors,
tors_from,
tors_to,
for_print=for_print)
return fig
[docs] def create_torsions_plot(self,
fig,
all_tors,
tors_from,
tors_to,
for_print=False):
"""
Creates a plot for torsions using the given matplot figure.
:param fig: Figure to draw the torsion subplots in.
:type fig: `matplotlib.figure.Figure`
:param tors: a List of List of torsions from each sublimation leg.
:type tors: list[list[[fep_edge_data.FEPTorsions]]
:param tors_from: Starting torsion index.
:type tors_from: int
:param tors_to: Ending torsion index
:type tors_to: int
:param for_print: Whether the figure is used for print.
:type for_print: bool
:returns: A list of the energy line plot axes
:rtype: list[matplotlib.axes.Axes]
"""
fig.clear()
nsub = len(all_tors)
gs = gridspec.GridSpec(tors_to - tors_from,
nsub,
top=0.90,
bottom=0.1,
left=0.05,
right=0.95,
wspace=0.1,
hspace=0.3) #, width_ratios=[3., 3.])
torsion_colors = []
line_plots = []
for irow in range(tors_from, tors_to):
# show tick labels location for the last row of torsion potential
if irow == tors_from:
hist_tick_pos = 'top'
elif irow == tors_to - 1:
hist_tick_pos = 'bottom'
else:
hist_tick_pos = None
for icol, leg_name in zip(range(nsub), self.fed.leg2_names_short):
tor = all_tors[icol].all_tors[irow]
tor_ax = fig.add_subplot(gs[irow % self.TORSIONS_PER_PAGE,
icol])
color = rh.colors[irow % len(rh.colors)]
torsion_colors.append(QtGui.QColor(color))
if icol == 0:
pot_tick_pos = 'left'
elif icol == nsub - 1:
pot_tick_pos = 'right'
else:
pot_tick_pos = None
tor_line_plot = tor.plot(tor_ax,
color=color,
for_print=for_print,
pot_tick_pos=pot_tick_pos,
hist_tick_pos=hist_tick_pos)
# set sublimatin leg name for each column
if hist_tick_pos == 'top':
tor_line_plot.title.set_y(1.5)
tor_line_plot.title.set_text(leg_name)
line_plots.append(tor_line_plot)
return line_plots
def _getTorsionStructures(self):
"""
:return: the two structures for the torsion plots
"""
return self.fed.mol_st, self.fed.mol_st
def _report_fep_simulation_details(self):
Ele = self.Elements
fed = self.fed
rh.pargph(Ele, '<u>System Details</u>')
rh.add_spacer(Ele)
rh.pargph(Ele, f"<font color='#888888'>Job name</font>: {fed.jobname}")
rh.pargph(
Ele,
f"<font color='#888888'>Job type</font>: {fed.perturbation_type}")
rh.pargph(
Ele, f"<font color='#888888'>Molecule Name:</font> {fed.mol_name}")
rh.add_spacer(Ele)
hyd_temp_str = '%0.1f' % fed.hyd_temperature
sub_temp_str = '%0.1f' % fed.sub_temperature
hyd_atoms_water = f'{fed.hyd_total_atoms} / {fed.hyd_total_waters}'
sub_atoms_water = f'{fed.sub_total_atoms} / {fed.sub_total_waters}'
hyd_delta_g = rh.get_pargph("%.2f ± %.2f" %
(fed.hyd_delta_g.val, fed.hyd_delta_g.unc),
fontsize=10)
delta_g_text = rh.get_pargph("ΔG [kcal/mol]",
fontsize=10,
color='gray')
lambda_text = rh.get_pargph("No. λ-Win.",
fontsize=10,
color='gray')
table = [[
'', 'Ensemble', 'Force Field', 'Temp. [K]', 'Sim. Time [ns]',
lambda_text, 'No. Atoms/Wats.', 'ASL Mol.', delta_g_text
]]
table.append([
f'{fed.leg1_name}', fed.hyd_ensemble, fed.hyd_ff, hyd_temp_str,
fed.hyd_sim_time, fed.hyd_total_replicas, hyd_atoms_water,
fed.hyd_mol_number, hyd_delta_g
])
for isub, sub in enumerate(fed.sub_delta_g):
sub_delta_g = rh.get_pargph("%.2f ± %.2f" % (sub[0], sub[1]),
fontsize=10)
table.append([
f'{fed.leg2_name}-{isub}', fed.sub_ensemble, fed.sub_ff,
sub_temp_str, fed.sub_sim_time, fed.sub_total_replicas,
sub_atoms_water, fed.sub_mol_number[isub], sub_delta_g
])
sub_delta_g = rh.get_pargph(
"%.2f ± %.2f" %
(fed.sub_median_delta_g.val, fed.sub_median_delta_g.unc),
fontsize=10)
table.append(
[f'{fed.leg2_name}-Med', '', '', '', '', '', '', '', sub_delta_g])
width = [1.2, 0.7, 0.7, 0.7, 0.9, 0.85, 1.15, 0.6, 1.1]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (1, 0), (-1, -1), 'CENTER'),
('ALIGN', (0, 1), (0, nrows - 1), 'RIGHT'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray),
('TEXTCOLOR', (0, 0), (0, nrows - 1), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.header(
Ele, "Transfer Free Energy (ΔΔG) "
"is: %.2f ± %.2f kcal/mol" %
(fed.solubility_dg.val, fed.solubility_dg.unc))
def _molecule_perturbation_details(self):
Ele = self.Elements
fed = self.fed
rh.pargph(Ele, '<u>Molecule Information</u>')
table = [[
'Name', 'HexID', 'No. Atoms / Heavy', 'Rot. Bonds', 'Atomic Mass',
'Charge', 'Mol. Formula'
]]
table.append([
fed.mol_name, fed.mol_hash,
f'{fed.mol_total_atoms} / {fed.mol_total_heavy}',
fed.mol_total_rot_bonds, f'{fed.mol_atomic_mass:.1f} au',
f'{"+" if fed.mol_charge > 0 else ""}{fed.mol_charge}',
fed.mol_formula
])
width = [2.2, 0.7, 1.1, 1.0, 1.0, 0.6, 1.2]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.add_spacer(Ele)
mol_2d_image = self.get_2d_ligand_image(fed.st)
table2 = [['Molecule'], [mol_2d_image]]
width2 = [4.0]
nfields2 = len(table2[0])
# yapf: disable
style2 = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields2 - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table2, style2, width2)
def _molecule_environment_report(self):
"""
Generate reports about each of the sublimation leg's alchemical
molecule and its environment
"""
Ele = self.Elements
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Molecular Environment of Sublimation Legs")
img = self._get_mol_env_plot(for_print=True)
Ele.append(img)
rh.pargph(Ele, self._get_mol_env_text())
def _get_mol_torsion_text(self):
text = """\
Rotatable bonds (Rb) in the molecule are enumerated and color-coded.
For each Rb, a representative dihedral angle is monitored throughout the
hydration and sublimation legs, its distribution is then plotted.
Hollow bars show <b>hydration</b> and filled bars show
<b>sublimation</b> leg distributions. Input starting conformation is
marked as a <font color="gray">gray</font> vertical line.
Potential energy around each Rb overlays the plot with the dark-blue
curve and corresponding labels on the Y-axis. Local strain energies are
shown below the plot. The units are in <i>kcal/mol</i>.
"""
return text
def _get_mol_env_text(self):
return """
For every sublimation leg, the interactions that the alchemical molecule
makes with its environment are monitored. These interactions are tracked
only for the replica window where the alchemical molecule is fully
coupled to the system. The interactions that are monitored are:
hydrogen bonds; intramolecular hydrogen (iH) bonds; halogen bonds;
π-π; π-cation and hydrophobic interactions. <i>Note:</i> to
improve the clarity and the readability of the plot, the number of
hydrophobic interactions are scaled down by 10x.
"""
def _get_mol_env_plot(self, for_print=True):
"""
:param for_print: Bool if the figure will be used for print
:type for_print: bool
"""
fig = figure.Figure(figsize=(6.0, 5.0))
FigureCanvas(fig)
self.create_mol_env_fig(fig, for_print=for_print)
temp_plot = self._get_temp_image_fn()
fig.savefig(temp_plot, bbox_inches='tight', dpi=300)
size = fig.get_size_inches()
(sx, sy) = rh.aspectScale(size[0], size[1], 6.0, 5.0)
return rh.platypus.Image(temp_plot, sx * rh.inch, sy * rh.inch)
[docs] def create_mol_env_fig(self, fig, for_print=False):
fontsize = 8 if for_print else 10
fig.clear()
ax = fig.add_subplot(111)
data = self.fed.sub_interactions(stats=True)
labels = self.fed.leg2_names_short
names = [
'Cation-Pi', 'H-Bond', 'Hal-Bond', 'Hydrphb. Inter.', 'Pi-Pi',
'Ionic Bond', 'iH-Bond'
]
keys = [
'CatPiResult', 'HBondResult', 'HalBondResult', 'HydrophobResult',
'PiPiResult', 'PolarResult', 'IntraHBondResult'
]
bottom = numpy.zeros(len(labels))
colors = reversed(rh.colors)
for key, name, color in zip(keys, names, colors):
vals = numpy.array([v[key][0] for v in data])
if not numpy.any(vals):
continue
if key == 'HydrophobResult':
vals /= 10.0
ax.bar(labels,
vals,
0.5,
label=name,
bottom=bottom,
color=color,
align='center')
bottom += vals
ax.set_ylabel('Number of Interactions', size=fontsize)
ax.yaxis.set_major_locator(MaxNLocator(nbins=5, prune='lower'))
_, ymax = ax.get_ylim()
# round up to the nearest 10s
ymax = math.ceil((ymax / 10.0)) * 10
ax.set_ylim((0, ymax))
ax.legend(bbox_to_anchor=(0., 1.02, 1.0, 0.102),
loc='lower left',
ncol=4,
mode="expand",
borderaxespad=0.,
title='Interaction Types',
fancybox=True,
prop={'size': 10})
if for_print:
rh.change_plot_colors(ax)
def _molecule_report(self):
"""
Generate ligand report
"""
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Compound Analysis for Physical λ-Replica")
table = [[
'Title ', 'PDB Name', 'No. Atoms', 'No. Heavy', 'Rot. Bonds',
'Atomic Mass', 'Charge'
]]
table.append([
fed.mol_name, fed.mol_pdb_name, fed.mol_total_atoms,
fed.mol_total_heavy, fed.mol_total_rot_bonds,
f'{fed.mol_atomic_mass:.1f} au',
f'{"+" if fed.mol_charge > 0 else ""}{fed.mol_charge}'
])
width = [1.5] + 6 * [0.85]
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER'),
('TEXTCOLOR', (0, 0), (nfields - 1, 0), rh.gray)]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.add_spacer(Ele)
self._mols_rmsd()
self._mols_sasa()
self._mols_molsa()
self._mols_psa()
self._mol_properties_text()
def _mol_properties_text(self):
text = """ \
Several molecular properties are tracked for each alchemical molecule
in all solubility legs (One hydration and numerous sublimation legs).
Here we show the time series and statistics for these properties. The
RMSD values are measured with respect to the input conformation.
"""
rh.pargph(self.Elements, text, fontsize=10)
def _mol_properties(self, header: str, data, axis_label):
"""
"""
Ele = self.Elements
rh.pargph(Ele, f'<u>{header}</u>')
rh.add_spacer(Ele)
psa_img = self.get_mol_series_plot(data, axis_label)
Ele.append(psa_img)
def _mols_psa(self):
fed = self.fed
header = 'Polar Surface Area'
numb_sub_legs = len(fed.sub_sid_psa(stats=False))
times = [fed.hyd_sid_snapshot_times_ns
] + [fed.sub_sid_snapshot_times_ns] * numb_sub_legs
psa = [fed.hyd_sid_psa(stats=False)] + fed.sub_sid_psa(stats=False)
names = [fed.leg1_name[:4]] + fed.leg2_names_short
data = list(zip(times, psa, names))
axis_label = r"PolSA ($\AA^2$)"
self._mol_properties(header, data, axis_label)
def _mols_molsa(self):
fed = self.fed
header = 'Molecular Surface Area'
numb_sub_legs = len(fed.sub_sid_molsa(stats=False))
times = [fed.hyd_sid_snapshot_times_ns
] + [fed.sub_sid_snapshot_times_ns] * numb_sub_legs
molsa = [fed.hyd_sid_molsa(stats=False)
] + fed.sub_sid_molsa(stats=False)
names = [fed.leg1_name[:4]] + fed.leg2_names_short
data = list(zip(times, molsa, names))
axis_label = r"MolSA ($\AA^2$)"
self._mol_properties(header, data, axis_label)
def _mols_sasa(self):
fed = self.fed
header = 'Solvent Accessible Surface Area'
numb_sub_legs = len(fed.sub_sid_sasa(stats=False))
times = [fed.hyd_sid_snapshot_times_ns
] + [fed.sub_sid_snapshot_times_ns] * numb_sub_legs
sasa = [fed.hyd_sid_sasa(stats=False)] + fed.sub_sid_sasa(stats=False)
names = [fed.leg1_name[:4]] + fed.leg2_names_short
data = list(zip(times, sasa, names))
axis_label = r"SASA ($\AA^2$)"
self._mol_properties(header, data, axis_label)
def _mols_rmsd(self):
fed = self.fed
header = 'Root Mean Square Deviation'
numb_sub_legs = len(fed.sub_sid_rmsd(stats=False))
times = [fed.hyd_sid_snapshot_times_ns
] + [fed.sub_sid_snapshot_times_ns] * numb_sub_legs
rmsds = [fed.hyd_sid_rmsd(stats=False)] + fed.sub_sid_rmsd(stats=False)
names = [fed.leg1_name[:4]] + fed.leg2_names_short
data = list(zip(times, rmsds, names))
axis_label = r"Molecules RMSD ($\AA$)"
self._mol_properties(header, data, axis_label)
[docs] def get_mol_series_plot(self, data, ylabel, for_print=True):
"""
:param data: a list of tuples for scatter plot
[(x-vals, y-vals, lables)...]
:type data: List(Tuple)`
:param ylabel: Label that goes on the y-axis
:type ylabel: str
:param for_print: Bool if the figure will be used for print
:type for_print: bool
"""
fig = figure.Figure(figsize=(7.2, 1.5))
FigureCanvas(fig)
self.create_mol_series_fig(data, fig, ylabel, for_print=for_print)
temp_plot = self._get_temp_image_fn()
fig.savefig(temp_plot, bbox_inches='tight', dpi=300)
size = fig.get_size_inches()
(sx, sy) = rh.aspectScale(size[0], size[1], 7.8, 1.8)
return rh.platypus.Image(temp_plot, sx * rh.inch, sy * rh.inch)
[docs] def create_mol_series_fig(self, data, fig, ylabel, for_print=False):
fontsize = 8 if for_print else 10
fig.clear()
spec = gridspec.GridSpec(ncols=2, nrows=1, width_ratios=[3, 2])
ax1 = fig.add_subplot(spec[0])
ax2 = fig.add_subplot(spec[1])
fig.subplots_adjust(wspace=0.05)
colors = [
'#025bac', '#df622a', '#a87458', '#de8d39', '#e2b481', '#896b31',
'#a88c2d', '#e3c647'
]
labels = []
for (ts, vals, name), color in zip(data, colors):
if not vals:
continue
p, = ax1.plot(ts, vals, color=color, label=name, alpha=0.8)
labels.append(name)
if ax1.get_ylim()[1] < 10:
ax1.set_ylim(ymin=0)
ax1.set_xlabel("Time (nsec)", size=fontsize)
ax1.set_ylabel(ylabel, size=fontsize)
ax1.yaxis.set_major_locator(MaxNLocator(nbins=4, prune='lower'))
# calculate and plot statistics
mean = [numpy.mean(d[1]) for d in data]
std = [numpy.std(d[1]) for d in data]
nbars = len(mean)
ax2.bar(range(nbars),
mean,
yerr=std,
color=colors[:nbars],
align='center',
ecolor='gray',
capsize=2)
ax2.set_ylim(ax1.get_ylim())
ax2.set_xticks(range(nbars))
ax2.set_xticklabels(labels)
ax2.yaxis.tick_right()
ax2.yaxis.set_major_locator(MaxNLocator(nbins=4, prune='lower'))
if for_print:
rh.change_plot_colors(ax1)
rh.change_plot_colors(ax2)
def _free_energy_convergence_profiles(self):
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Free Energy Convergence")
rh.pargph(Ele, '<u>%s Leg</u>' % fed.leg1_name)
rh.add_spacer(Ele)
hyd_convergence = self._gen_free_energy_convergence_plot(
fed.hyd_start_time_ns, fed.hyd_end_time_ns, fed.hyd_delta_g_forward,
fed.hyd_delta_g_forward_err, fed.hyd_delta_g_reverse,
fed.hyd_delta_g_reverse_err, fed.hyd_delta_g_sliding,
fed.hyd_delta_g_sliding_err)
Ele.append(hyd_convergence)
rh.add_spacer(Ele)
# sublimaiton leg plots
for isub, (sub_delta_g_forward, sub_delta_g_forward_err,
sub_delta_g_reverse, sub_delta_g_reverse_err,
sub_delta_g_sliding, sub_delta_g_sliding_err) in enumerate(
zip(fed.sub_delta_g_forward, fed.sub_delta_g_forward_err,
fed.sub_delta_g_reverse, fed.sub_delta_g_reverse_err,
fed.sub_delta_g_sliding,
fed.sub_delta_g_sliding_err)):
# plot 3 plots per page
if not (isub + 1) % 3:
rh.pargph(Ele,
self._get_free_energy_convergence_text(table=False))
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Free Energy Convergence (Cont.)")
rh.pargph(Ele, f'<u>{fed.leg2_name}-{isub} Leg</u>')
rh.add_spacer(Ele)
sub_convergence = self._gen_free_energy_convergence_plot(
fed.sub_start_time_ns, fed.sub_end_time_ns, sub_delta_g_forward,
sub_delta_g_forward_err, sub_delta_g_reverse,
sub_delta_g_reverse_err, sub_delta_g_sliding,
sub_delta_g_sliding_err)
Ele.append(sub_convergence)
rh.add_spacer(Ele)
rh.pargph(Ele, self._get_free_energy_convergence_text(table=False))
def _replica_exchange_density(self):
Ele = self.Elements
fed = self.fed
rh.new_page(Ele)
rh.report_add_logo(Ele)
rh.header(Ele, "Exchange Density of FEP Replicas Over "
"λ-Windows")
rh.pargph(Ele, f'<u>{fed.leg1_name} Leg</u>')
rh.add_spacer(Ele)
hyd_img = self.get_rest_density_img(fed.hyd_rest_exchanges)
Ele.append(hyd_img)
rh.add_spacer(Ele)
rh.pargph(Ele, f'<u>{fed.leg2_name} Legs</u>')
rh.add_spacer(Ele)
table = []
total_sub_legs = len(fed.sub_rest_exchanges)
plots_per_row = 3
for irow in range(math.ceil(total_sub_legs / plots_per_row)):
name_row = []
plot_row = []
for icol in range(3):
iplot = irow * plots_per_row + icol
if iplot < total_sub_legs:
name_row.append(f'{fed.leg2_name}-{iplot}')
plot_row.append(
self.get_rest_density_img(fed.sub_rest_exchanges[iplot],
legend=False,
size_ratio=0.53))
else:
name_row.append('')
plot_row.append('')
table.append(name_row)
table.append(plot_row)
width = [2.6] * 3
nfields = len(table[0])
nrows = len(table)
# yapf: disable
style = [('BOTTOMPADDING', (0, 1), (-1, -1), 1),
('TOPPADDING', (0, 1), (-1, -1), 1),
('ALIGN', (0, 0), (-1, -1), 'CENTER')]
# yapf: enable
rh.add_table(Ele, table, style, width)
rh.pargph(Ele, self.get_rest_density_text, fontsize=10)