"""
Functions for creating or converting pose viewer files.
The module can convert a 'pose viewer' type file into a series of complexes,
and convert complexes into ligand-only, receptor-only, or pose viewer files.
The script pv_convert.py is the preferred command-line interface to this module.
Copyright Schrodinger, LLC. All rights reserved.
"""
################################################################################
# Packages
################################################################################
import warnings
import schrodinger.structutils.analyze as analyze
import schrodinger.utils.log as log
from schrodinger import structure
from schrodinger.application.glide import MAX_LIGAND_ATOMS
################################################################################
# Logging
################################################################################
logger = log.get_output_logger(__file__)
logger.setLevel(log.logging.INFO)
################################################################################
# Functions
################################################################################
[docs]def is_pv_file(
file_name,
max_ligand_atoms=MAX_LIGAND_ATOMS,
test_all=True,
test_recep_atom_count=True,
test_ligand_atom_count=True,
test_recep_gscore=True,
test_ligand_gscore=True,
):
"""
Returns a True if the file appears to be a pose view format file.
The tests are simple, and may not be conclusive. It is assumed that
the trivial test for a 'pv.mae' extension has already been performed.
Specific tests can be skipped by setting the test_recep_atom_count,
test_recep_gscore, test_ligand_atom_count, and test_ligand_gscore
booleans.
:type max_ligand_atoms: int
:param max_ligand_atoms:
The maximum size in atoms of a putative ligand, and minimum
size of the receptor. The default is the package constant
schrodinger.application.glide.MAX_LIGAND_ATOMS.
:type test_all: bool
:param test_all:
Perform all tests, regardless of other test_* booleans.
Default = True.
:type test_recep_atom_count: bool
:param test_recep_atom_count:
Receptor must have more than max_ligand_atoms. Default =
True.
:type test_ligand_atom_count: bool
:param test_ligand_atom_count:
Ligand must have less than or equal to max_ligand_atoms.
Default = True.
:type test_recep_gscore: bool
:param test_recep_gscore:
Receptor must not have 'r_i_glide_gscore' property.
Default = True.
:type test_ligand_gscore: bool
:param test_ligand_gscore:
Ligand must have 'r_i_glide_gscore' property. Default = True.
"""
msg = ("is_pv_file is a deprecated function; use is_valid_pv_file in"
"schrodinger.application.glide.utils instead")
warnings.warn(msg, DeprecationWarning, stacklevel=2)
putative_receptor = None
putative_poses = []
# Fetch some structures to assay
logger.info("Sampling structures to assay pv format.")
for index, st in enumerate(structure.StructureReader(file_name)):
if index == 0:
putative_receptor = st
if index in [1, 2, 3, 5, 10]:
putative_poses.append(st)
if index > 10:
break
# Receptor should have more than 'max_ligand_atoms'
if test_all or test_recep_atom_count:
if not len(putative_receptor.atom) > max_ligand_atoms:
logger.info("Receptor has too few atoms.")
return False
# Ligand should have fewer that 'max_ligand_atoms'
if test_all or test_ligand_atom_count:
for st in putative_poses:
if len(st.atom) > max_ligand_atoms:
logger.info("Putative ligand pose has too many atoms.")
return False
# Receptor should not have a glide score. Optional test because
# complexes may be merged and split such that the 'receptor' contains
# pose properties, but it is not likely.
if test_all or test_recep_gscore:
if putative_receptor.property.get('r_i_glide_gscore'):
logger.info("Putative receptor has a glide score.")
return False
# Ligands should have glide score, this is an optional test because
# you may need pv format, but not require glide props. e.g. MBAE,
# or prime_mmgbsa input
if test_all or test_ligand_gscore:
for st in putative_poses:
if not st.property.get('r_i_glide_gscore'):
logger.info("Putative ligand pose is missing glide score.")
return False
# File has passed all tests, it could be pv format
return True
##########################################################################
# Classes
##########################################################################
[docs]class Complex:
"""
A helper class to split a receptor-ligand complex structure into
ligand, and receptor structure components, and write the structure
components to new files.
API Examples::
# Write a pv file
comp = Complex(st)
comp.writePv('output_pv.mae')
# Write a ligand file, transfer all the complex's properties
comp = Complex(
st
ligand_properties = st.property.keys()
)
comp.writeLigand('ligand-only.mae')
The receptor will have the same properties as the input complex.
By default, the ligand takes the title of the complex but has no
other properties. The ligand title can be selected from one of
the complex's properties with ligand_title_source, and a list
of properties to transfer to the ligand can be specified with
ligand_properties.
:vartype complex: structure.Structure
:ivar complex:
The original receptor-ligand complex.
:vartype receptor: structure.Structure
:ivar receptor:
The receptor-only (non-ligand).
:vartype ligand: structure.Structure
:ivar ligand:
The ligand-only structure.Structure.
:vartype ligand_indexes: list
:ivar ligand_indexes:
List of atom numbers for the ligand in the complex, the
match to the ASL evaluation.
:vartype ligand_title_source: string
:ivar ligand_title_source:
A single m2io dataname of an entry level property in the
complex structure to use as the title for the ligand.
Default is 's_m_title'.
:vartype ligand_properties: list
:ivar ligand_properties:
A list of m2io datanames of entry level properties to copy
from the complex to the ligand structure.
:raise RuntimeError: If the structure can't be split.
"""
[docs] def __init__(self,
st,
ligand_asl=None,
ligand_title_source=None,
ligand_properties=None):
"""
:type st: structure.Structure
:param st:
Structure instance of the receptor-ligand complex.
:type ligand_asl: str
:param ligand_asl:
Optional Atom Selection Language expression for the ligand
molecule, Default is None, which will select the last molecule
in the complex (mol.n #).
:type ligand_title_source: str
:param ligand_title_source:
An m2io entry data name of the complex property to use as the
ligand's title. Default is 's_m_title'.
:type ligand_properties: list of str
:param ligand_properties:
List of m2io entry level data names in the complex to be
transfered to the ligand entry.
"""
self.complex = st
self.receptor = None
self.ligand = None
self.ligand_indexes = None
self.ligand_asl = None
# Check to see we have a split-worthy molecule
if not ligand_asl and self.complex.mol_total < 2:
raise RuntimeError(
'Unable to detect ligand: no ligand ASL'
' specified and structure has only one logical molecule.')
if ligand_asl:
self.ligand_asl = ligand_asl
else:
self.ligand_asl = 'mol.n %d' % self.complex.mol_total
if ligand_title_source:
self.ligand_title_source = ligand_title_source
else:
self.ligand_title_source = 's_m_title'
if ligand_properties:
self.ligand_properties = ligand_properties
else:
self.ligand_properties = []
# Do the split
self._split()
# FIXME This is a good, widely applicable asl, consider if/how it
# should be incorporated.
# self.ligand_asl = '(mol.atoms 5-130) AND NOT (( backbone or atom.pt 2HA, 3HA ) OR ( sidechain ) )'
return
def _split(self):
"""
Extracts the ligand from the ligand-receptor complex, creates
the ligand structure attribute, and deletes the ligand from the
complex to form the receptor attribute.
"""
# Evaluate ASL, make sure the match is unambiguous
self.ligand_indexes = analyze.evaluate_asl(self.complex,
self.ligand_asl)
mol_set = set()
for index in self.ligand_indexes:
mol_set.add(self.complex.atom[index].molecule_number)
if len(mol_set) != 1:
msg = 'ASL "{0}" identifies'.format(self.ligand_asl)
if len(mol_set) > 1:
msg += ' atoms from more than one molecule.'
else:
msg += ' no atoms.'
raise RuntimeError(msg)
# Extract the ligand, set properties
self.ligand = self.complex.extract(self.ligand_indexes)
for prop in self.ligand_properties:
if prop in self.complex.property:
self.ligand.property[prop] = self.complex.property[prop]
# Title as separate property, and done second to avoid clobbering
self.ligand.title = self.complex.property[self.ligand_title_source]
# Create the receptor structure
self.receptor = self.complex.copy()
self.receptor.deleteAtoms(self.ligand_indexes)
return
[docs] def writePv(self, file_name):
"""
Writes receptor and ligand to disk, in pv file format.
"""
self.receptor.write(file_name)
self.ligand.append(file_name)
return
[docs] def writeLigand(self, file_name):
"""
Writes the ligand structure to the passed filename, clobbering
as needed.
"""
self.ligand.write(file_name)
return
[docs] def appendLigand(self, file_name):
"""
Appends the ligand structure to the passed filename.
"""
self.ligand.append(file_name)
return
[docs] def writeReceptor(self, file_name):
"""
Writes the receptor structure to the passed filename, clobbering
as needed.
"""
self.receptor.write(file_name)
return
[docs] def appendReceptor(self, file_name):
"""
Appends the receptor structure to the passed filename.
"""
self.receptor.append(file_name)
return
# EOF