import schrodinger.structutils.analyze as analyze
import schrodinger.structutils.build as build
import schrodinger.structutils.smiles as smiles
import schrodinger.utils.sea as sea
from schrodinger.application.bioluminate import protein
from schrodinger.application.desmond import constants
from schrodinger.application.desmond.struc import get_reference_ct
PROTEIN_ASL = "(protein)"
MAX_RING_SIZE = 8
[docs]def removeAminoAcids(ligand_list):
"""
This function iterates over ligand structures in the given list
and removes those that have amino acid names in their pdb residue
name.
:param ligand_list: list of ligand structures
:type ligand_list: list
:return: a list of `Ligand` instances
:rtype: list
"""
AMINO_ACID_NAMES = (protein.Mutator.SUPPORTED_BUILD_RESIDUES +
["HID", "ARN", "ASH", "GLH", "LYN"])
new_list = []
for ligand in ligand_list:
if not ligand.st.atom[1].pdbres.strip() in AMINO_ACID_NAMES:
new_list.append(ligand)
return new_list
[docs]def getLigand(cms_st):
"""
This parses a CMS for the ligand to use.
:param cms_st: A CMS to find a ligand within
:type cms_st: cms.Cms object
:return: (ligand structure, ligand asl)
:rtype: (Structure, str)
"""
ligand_list = analyze.find_ligands(cms_st)
# only filter ligand list when it has more than 1 ligand
if len(ligand_list) > 1:
ligand_list = removeAminoAcids(ligand_list)
if len(ligand_list) == 0:
return None, None
ligand_st = ligand_list[0].st
ligand_asl = str(ligand_list[0].ligand_asl)
return ligand_st, ligand_asl
[docs]def getASLExcludingLigand(asl_type, protein_asl=PROTEIN_ASL, ligand_asl=''):
"""
This gives an ASL for a subsection of the protein without the ligand.
:param asl_type: Type of protein wanted, Heavy, Backbone, etc
:type asl_type: str
:param protein_asl: The ASL to describe the entire protein
:type protein_asl: str
:param ligand_asl: The ASL to describe the ligand
:type ligand_asl: str
:return: ASL for the protein component minus the ligand
:rtype: str
"""
if ligand_asl:
exclude_ligand = 'and not (%s)' % ligand_asl
else:
exclude_ligand = ''
types = {
'Heavy': "(%s and not (atom.ele H) %s)",
'Backbone': "((%s and backbone) and not (atom.ele H) %s)",
'C-Alpha': "((%s and backbone ) and ( atom.ptype ' CA ') %s)",
'All residues': "(%s %s)",
'Side chains': "(%s and sidechain %s)"
}
return types[asl_type] % (protein_asl, exclude_ligand)
[docs]def getRMSDKeywords(protein_asl, ligand_asl, ref_struct_fname=None, frame=0):
"""
Returns keywords for RMSD analysis of protein components.
:param protein_asl: The ASL to describe the entire protein
:type protein_asl: str
:param ligand_asl: The ASL to describe the ligand
:type ligand_asl: str or None
:param ref_struct_fname: Path to the structure to do the RMSD against
:type ref_struct_fname: None or `str`
:param frame: The frame to take the RMSD against, ignored if
ref_struct_fname!=None
:type frame: int
:return: Keyword(s) for calculation
:rtype: sea.Map
"""
types = ["Backbone", "Side chains", "C-Alpha", "Heavy"]
for rmsd_type in types:
asl = getASLExcludingLigand(rmsd_type, protein_asl, ligand_asl)
ark_kw = sea.Map()
ark_kw["RMSD"] = sea.Map()
ark_kw["RMSD"]["Unit"] = "Angstrom"
#Either use a reference frame, or a reference structure
if ref_struct_fname:
ark_kw["RMSD"]["Reference_Struct"] = ref_struct_fname
else:
ark_kw["RMSD"]["Frame"] = frame
ark_kw["RMSD"]["Type"] = "ASL"
ark_kw["RMSD"]["SelectionType"] = rmsd_type
ark_kw["RMSD"]["ASL"] = asl
ark_kw["RMSD"]["Panel"] = 'pl_interact_survey'
ark_kw["RMSD"]["Tab"] = 'pl_rmsd_tab'
if asl:
ark_kw["RMSD"]["ASL"] = asl
yield ark_kw
[docs]def getRMSDLigandKW(ligand_asl, fitby, ref_struct_fname=None, frame=0):
"""
Returns keywords for RMSD analysis of the ligand.
:param ligand_asl: The ASL to describe the ligand
:type ligand_asl: str
:param fitby: None or an ASL to describe what portion of the CMS to fit
:type fitby: None or str
:param ref_struct_fname: Path to the structure to do the RMSD against
:type ref_struct_fname: None or `str`
:param frame: The frame to take the RMSD against, ignored if
ref_struct_fname!=None
:type frame: int
"""
ark_kw = sea.Map()
ark_kw["RMSD"] = sea.Map()
ark_kw["RMSD"]["Unit"] = "Angstrom"
#Either use a reference frame, or a reference structure
if ref_struct_fname:
ark_kw["RMSD"]["Reference_Struct"] = ref_struct_fname
else:
ark_kw["RMSD"]["Frame"] = frame
ark_kw["RMSD"]["Type"] = "Ligand"
ark_kw["RMSD"]["SelectionType"] = "Ligand"
ark_kw["RMSD"]["Panel"] = 'pl_interact_survey'
ark_kw["RMSD"]["Tab"] = 'pl_rmsd_tab'
ark_kw["RMSD"]["UseSymmetry"] = 'yes'
ark_kw["RMSD"]["ASL"] = ligand_asl
if fitby:
ark_kw["RMSD"]["FitBy"] = fitby
return ark_kw
[docs]def getRMSFProtKeywords(protein_asl,
ligand_asl,
ref_struct_fname=None,
frame=0):
"""
Returns keywords for RMSF analysis of protein components.
:param protein_asl: The ASL to describe the entire protein
:type protein_asl: str
:param ligand_asl: The ASL to describe the ligand
:type ligand_asl: str
:param ref_struct_fname: A structure to do the RMSF against
:type ref_struct_fname: None or schrodinger.structure.Structure
:param frame: The frame to take the RMSF against, ignored if
ref_struct_fname!=None
:type frame: int
:return: Keyword(s) for calculation
:rtype: sea.Map
"""
types = ["Backbone", "Side chains", "C-Alpha", "Heavy"]
for rmsf_type in types:
asl = getASLExcludingLigand(rmsf_type, protein_asl, ligand_asl)
ark_kw = sea.Map()
ark_kw["RMSF"] = sea.Map()
ark_kw["RMSF"]["Unit"] = "Angstrom"
#Either use a reference frame, or a reference structure
if ref_struct_fname:
ark_kw["RMSF"]["Reference_Struct"] = ref_struct_fname
else:
ark_kw["RMSF"]["Frame"] = frame
#ark_kw["RMSF"]["FitBy"] = "backbone"
ark_kw["RMSF"]["FitBy"] = getASLExcludingLigand("Backbone", protein_asl,
ligand_asl)
ark_kw["RMSF"]["Type"] = "ASL"
ark_kw["RMSF"]["SelectionType"] = rmsf_type
ark_kw["RMSF"]["ASL"] = asl
ark_kw["RMSF"]["Panel"] = 'pl_interact_survey'
ark_kw["RMSF"]["Tab"] = 'p_rmsf_tab'
yield ark_kw
[docs]def getSSEProtKeywords(protein_asl, ligand_asl):
"""
Returns keywords for SSE analysis
:param protein_asl: The ASL to describe the entire protein
:type protein_asl: str
:param ligand_asl: The ASL to describe the ligand
:type ligand_asl: str
:return: Keyword(s) for calculation
:rtype: sea.Map
"""
# add secondary structure analysis for protein
ark_kw = sea.Map()
ark_kw["SecondaryStructure"] = sea.Map()
asl = getASLExcludingLigand("All residues", protein_asl, ligand_asl)
ark_kw["SecondaryStructure"]["Type"] = "ASL"
ark_kw["SecondaryStructure"]["SelectionType"] = asl
ark_kw["SecondaryStructure"]["ASL"] = asl
ark_kw["SecondaryStructure"]["Panel"] = "pl_interact_survey"
ark_kw["SecondaryStructure"]["Tab"] = "p_sse_tab"
ark_kw["SecondaryStructure"][
"Legend"] = "-1: NONE; 0: LOOP; 1: HELIX; 2: STRAND; 3:TURN"
yield ark_kw
[docs]def getRMSFLigandKW(fitby_asl, ligand_asl, ref_struct_fname=None, frame=0):
"""
Returns keywords for RMSF analysis of the ligand.
:param ligand_asl: The ASL to describe the ligand
:type ligand_asl: str
:param fitby: None or an ASL to describe what portion of the CMS to fit
:type fitby: None or str
:param ref_struct_fname: A structure to do the RMSF against
:type ref_struct_fname: None or schrodinger.structure.Structure
:param frame: The frame to take the RMSF against, ignored if
ref_struct_fname!=None
:type frame: int
"""
ark_kw = sea.Map()
ark_kw["RMSF"] = sea.Map()
ark_kw["RMSF"]["Unit"] = "Angstrom"
#Either use a reference frame, or a reference structure
if ref_struct_fname:
ark_kw["RMSF"]["Reference_Struct"] = ref_struct_fname
else:
ark_kw["RMSF"]["Frame"] = frame
ark_kw["RMSF"]["FitBy"] = fitby_asl
ark_kw["RMSF"]["ASL"] = ligand_asl
ark_kw["RMSF"]["Type"] = "ASL"
ark_kw["RMSF"]["Panel"] = 'pl_interact_survey'
ark_kw["RMSF"]["Tab"] = 'l_rmsf_tab'
return ark_kw
[docs]def getProtLigInterKW(protein_asl, ligand_asl, metal_asl=None):
"""
Returns keywords for PLI analysis
:param protein_asl: The ASL to describe the entire protein
:type protein_asl: str
:param ligand_asl: The ASL to describe the ligand
:type ligand_asl: str
:param metal_asl: The ASL to describe metals or ions.
If None, use all metals and ions in
the input structure.
:type metal_asl: str or None
"""
ark_kw = sea.Map()
ark_kw["ProtLigInter"] = sea.Map()
ark_kw["ProtLigInter"]["ASL"] = protein_asl
ark_kw["ProtLigInter"]["LigandASL"] = ligand_asl
if metal_asl is not None:
ark_kw["ProtLigInter"]["MetalASL"] = metal_asl
ark_kw["ProtLigInter"]["Panel"] = 'pl_interact_survey'
ark_kw["ProtLigInter"]["Tab"] = 'pl_inter_tab'
return ark_kw
[docs]def getSolubilityKeywords(molecule_asl):
"""
Returns keywords for analyzing solubility's sublimation leg.
Here we look at the environment of the molecule and its interactions
with it.
"""
ark_kw = sea.Map()
ark_kw['AmorphousCrystalInter'] = sea.Map()
ark_kw['AmorphousCrystalInter']['ASL'] = molecule_asl
return ark_kw
[docs]def get_idx(ligand_st, lig_idx):
"""
This will get the atoms index associated with an asl index
:param ligand_st : ligand structure
:type ligand_st : structure
:param lig_idx : ligand atom index
:type lig_idx : int
:return : atom_index of the asl selection for that atom
:rtype : atom_index
"""
if constants.ORIGINAL_INDEX in list(ligand_st.atom[lig_idx].property):
atom_index = ligand_st.atom[lig_idx].property[constants.ORIGINAL_INDEX]
elif 'i_i_internal_atom_index' in list(ligand_st.atom[lig_idx].property):
atom_index = ligand_st.atom[lig_idx].property['i_i_internal_atom_index']
#FIXME: what if both condition is False?
return atom_index
[docs]def sort_atoms(st, a_id, exclude_atom_id=None, for_fep=False):
"""
Extract atoms that are bonded to a_id atom (excluding the exclude_atom_id
atom) and returns a list of atoms in the order such that hydrogens are last
:param st: small molecule structure
:type st: structure
:param a_id: index of the atom whose bonded atoms you want to return
:type a_id: int
:param exclude_atom_id: remove that atom from the list
:type exclude_atom_id: int
:param for_fep: if this options is ture, then return only atoms that is
mapped to another atoms (contains i_fep_mapping prop)
:type for_fep: bool
:rtype: int
:return: The index of the first heavy atom a_id, that is bounded, and that
doesn't exclude_atom_id
"""
bonded_atoms = []
# sort atoms such that hydrogen and/or linear atoms appear last
for atom in st.atom[a_id].bonded_atoms:
if atom.index == exclude_atom_id:
continue
a_type = 1 if atom.atomic_number == 1 else 0
bond_order = st.getBond(st.atom[a_id], atom).order
if bond_order == 0:
continue
cm = sum([a.atomic_weight for a in atom.bonded_atoms])
angle = st.measure(exclude_atom_id, a_id,
atom.index) if exclude_atom_id else 0.0
# this is to make the algorithm to be less conformation dependant but
# we want to bias the selection away from linear atoms.
angle = angle if angle > 170. else 0.
bonded_atoms.append((a_type, bond_order, cm, angle, atom.index))
bonded_atoms = [atom_props[-1] for atom_props in sorted(bonded_atoms)]
if not for_fep:
return bonded_atoms[0]
# if for_fep, make sure to return the mapped atom
for ba in bonded_atoms:
if st.atom[ba].property.get(constants.FEP_MAPPING, 0):
return ba
return None
[docs]def get_rotatable_bonds(st):
"""
returns all rotatable bonds, defined as torsions. returns original atoms indeces.
:param st: structure of a ligand
:type st: structure
:rtype: list
:return: a list of four atoms that define a rotatable bond
"""
rb_atoms = []
for rb in analyze.rotatable_bonds_iterator(st, max_size=MAX_RING_SIZE):
a2, a3 = rb
a1 = sort_atoms(st, a2, a3)
a4 = sort_atoms(st, a3, a2)
b1 = st.atom[a1].property[constants.ORIGINAL_INDEX]
b2 = st.atom[a2].property[constants.ORIGINAL_INDEX]
b3 = st.atom[a3].property[constants.ORIGINAL_INDEX]
b4 = st.atom[a4].property[constants.ORIGINAL_INDEX]
rb_atoms.append([b1, b2, b3, b4])
return rb_atoms
[docs]def canonicalize_ligand(st):
"""
Use uSMILES to canonicalize the ligand structure, this is so we get the
order of the rotatable bonds(torsions) in the same order, regardless of
the atom order
:param st: structure of a ligand
:type st: structure
:rtype: structure
:return: a structure file with reordered atoms according to unique SMILES
"""
smiles_generator = smiles.SmilesGenerator(
smiles.STEREO_FROM_ANNOTATION_AND_GEOM, unique=True)
# get all heavy atoms and polar hydrogens
asl = '(not a.e H) or ( atom.ele H and not /C0-H0/ )'
heavy_polarH_list = analyze.evaluate_asl(st, asl)
new_st = st.extract(heavy_polarH_list)
# get a list of left over hydrogens
pol_h = [atom.index for atom in new_st.atom if atom.element == 'H']
smi, atom_list = smiles_generator.getSmilesAndMap(new_st)
# add the polar hydrogens to the end of the list
for h in pol_h:
if h not in atom_list:
atom_list.append(h)
return build.reorder_atoms(new_st, atom_list)
[docs]def getLigandPropsKeywords(ligand_asl,
calcRMSD,
ref_struct_fname=None,
frame=None):
"""
Returns keywords for ligand-specific analysis
* Intramolecular Hbonds
* Molecular Surface Area (Cannolly surface)
* Solvent Accessible Surface Area
* Polar Surface area
* Radius of Gyration
:param ligand_asl: The ASL to describe the ligand
:type ligand_asl: str
:param calcRMSD: Bool to also set up RMSD jobs
:type calcRMSD: bool
:param ref_struct_fname: path to the reference structure
:param frame: frame index
:type frame: `int`
:return: Keyword(s) for calculation
:rtype: sea.Map
"""
ark_kw = sea.Map()
ark_kw["LigandHBonds"] = sea.Map()
ark_kw["LigandHBonds"]["ASL1"] = ligand_asl
ark_kw["LigandHBonds"]["Name"] = 'Intramolecular Hydrogen Bonds'
ark_kw["LigandHBonds"]["Type"] = 'ASL'
ark_kw["LigandHBonds"]["Unit"] = 'Numb. of HBonds'
ark_kw["LigandHBonds"]["Panel"] = 'pl_interact_survey'
ark_kw["LigandHBonds"]["ReturnHBonds"] = 'True'
ark_kw["LigandHBonds"]["Tab"] = 'l_properties_tab'
yield ark_kw
ark_kw = sea.Map()
ark_kw["Molecular_Surface_Area"] = sea.Map()
ark_kw["Molecular_Surface_Area"]['ASL'] = ligand_asl
ark_kw["Molecular_Surface_Area"]['Grid_Spacing'] = 0.5
ark_kw["Molecular_Surface_Area"]["Unit"] = 'Angstrom^2'
ark_kw["Molecular_Surface_Area"]['Panel'] = 'pl_interact_survey'
ark_kw["Molecular_Surface_Area"]['Tab'] = 'l_properties_tab'
yield ark_kw
ark_kw = sea.Map()
ark_kw["SA_Surface_Area"] = sea.Map()
ark_kw["SA_Surface_Area"]['ASL'] = ligand_asl
# TODO: need a better definition to exclude FEP image
ions_asl = '(metals) or (metalloids)'
if 'protein' in ligand_asl:
exclude_asl = 'not ((%s) or (ions) or (water) or %s)' % (ligand_asl,
ions_asl)
else:
exclude_asl = 'not ((protein) or (%s) or (ions) or (water) or %s)' % \
(ligand_asl, ions_asl)
ark_kw["SA_Surface_Area"]['Exclude_ASL'] = exclude_asl
ark_kw["SA_Surface_Area"]['Resolution'] = 0.3
ark_kw["SA_Surface_Area"]["Unit"] = 'Angstrom^2'
ark_kw["SA_Surface_Area"]['Panel'] = 'pl_interact_survey'
ark_kw["SA_Surface_Area"]['Tab'] = 'l_properties_tab'
yield ark_kw
ark_kw = sea.Map()
ark_kw["Polar_Surface_Area"] = sea.Map()
ark_kw["Polar_Surface_Area"]['ASL'] = ligand_asl
ark_kw["Polar_Surface_Area"]['Resolution'] = 0.3
ark_kw["Polar_Surface_Area"]["Unit"] = 'Angstrom^2'
ark_kw["Polar_Surface_Area"]['Panel'] = 'pl_interact_survey'
ark_kw["Polar_Surface_Area"]['Tab'] = 'l_properties_tab'
yield ark_kw
ark_kw = sea.Map()
ark_kw["Rad_Gyration"] = sea.Map()
ark_kw["Rad_Gyration"]['ASL'] = ligand_asl
ark_kw["Rad_Gyration"]["Unit"] = 'Angstrom'
ark_kw["Rad_Gyration"]['Panel'] = 'pl_interact_survey'
ark_kw["Rad_Gyration"]['Tab'] = 'l_properties_tab'
yield ark_kw
if calcRMSD:
ark_kw = sea.Map()
ark_kw["RMSD"] = sea.Map()
ark_kw["RMSD"]["Unit"] = "Angstrom"
if ref_struct_fname:
ark_kw["RMSD"]["Reference_Struct"] = ref_struct_fname
else:
ark_kw["RMSD"]["Frame"] = frame
ark_kw["RMSD"]["Type"] = "Ligand"
ark_kw["RMSD"]["SelectionType"] = "Ligand"
ark_kw["RMSD"]["Panel"] = 'pl_interact_survey'
ark_kw["RMSD"]["Tab"] = 'l_properties_tab'
ark_kw["RMSD"]["UseSymmetry"] = 'yes'
ark_kw["RMSD"]["ASL"] = ligand_asl
yield ark_kw
def _getLoneTorsionFEP(lig_st, core_rb):
"""
This method returns Torsions/RB that are not in the core_rb.
"""
ret_torsions = []
for rb in analyze.rotatable_bonds_iterator(lig_st, max_size=MAX_RING_SIZE):
a2, a3 = rb
if (a2, a3) in core_rb or (a3, a2) in core_rb:
continue
a1 = sort_atoms(lig_st, a2, a3)
a4 = sort_atoms(lig_st, a3, a2)
ret_torsions.append((a1, a2, a3, a4))
return ret_torsions
[docs]def getFEPTorsionKeywords(lig1_st,
lig2_st,
lig1_asl,
lig2_asl,
fep_lambda=0,
is_covalent=False):
"""
For a pair for ligands/fragments that, whose atoms are mapped, return sea
object of torsions for each rotatable bond.
"""
def get_ark_tor(lig_st, rb, asl, fep_pair=None):
fsys_ct_rb = lig_to_fullsys_aid(rb, lig_st)
# Use the analysis reference coords if present
ref_lig_st = get_reference_ct(lig_st)
ref_lig_st = get_reference_ct(
ref_lig_st,
prop_names=constants.ANALYSIS_REFERENCE_COORD_PROPERTIES)
ark_kw = sea.Map()
ark_kw["Torsion"] = sea.Map()
ark_kw["Torsion"]["a1"] = fsys_ct_rb[0]
ark_kw["Torsion"]["a2"] = fsys_ct_rb[1]
ark_kw["Torsion"]["a3"] = fsys_ct_rb[2]
ark_kw["Torsion"]["a4"] = fsys_ct_rb[3]
ark_kw["Torsion"]["Unit"] = "Degrees"
ark_kw["Torsion"]["Panel"] = 'pl_interact_survey'
ark_kw["Torsion"]["Tab"] = 'l_torsions_tab'
ark_kw["Torsion"]["LigandASL"] = asl
ark_kw["Torsion"]["StartingConformation"] = ref_lig_st.measure(*rb)
if fep_pair is not None:
ark_kw["Torsion"]["FEP_pair"] = fep_pair
return ark_kw
core1, core2, lone1, lone2 = _getFEPTorsion(lig1_st, lig2_st, is_covalent)
ret_list = sea.List()
if fep_lambda == 0:
for irb, rb in enumerate(core1, start=1):
ret_list.append(get_ark_tor(lig1_st, rb, lig1_asl, irb))
for rb in lone1:
ret_list.append(get_ark_tor(lig1_st, rb, lig1_asl))
else:
for irb, rb in enumerate(core2, start=1):
ret_list.append(get_ark_tor(lig2_st, rb, lig2_asl, irb))
for rb in lone2:
ret_list.append(get_ark_tor(lig2_st, rb, lig2_asl))
return ret_list
def _valid_mapped_fep_atom(atom):
"""
Check if atom is alchemically mapped
"""
return atom.property.get(constants.FEP_MAPPING, 0) != 0
def _getFEPTorsion(lig1_st, lig2_st, is_covalent):
"""
This method generates torsions that represent rotatable bonds by taking two
ligands, that are alchemically-mapped. The torsions generated are mapped
between the two ligands. CoreRBs are the ones where the dihedrals are
mapped, loneRBs are those that are not mapped.
:param is_covalent: if the covalent ligand is used, the provided structures
will be a subset -- sidechain residue + ligand
:type is_covalent: bool
"""
core1_rb = []
core2_rb = []
cpx1_frag1_map = None
if is_covalent:
cpx1_frag1_map = {
atom.property[constants.ORIGINAL_INDEX]: atom.index
for atom in lig1_st.atom
}
# Add mapping information to lig1_st
for atom in lig2_st.atom:
m_idx = atom.property.get(constants.FEP_MAPPING, 0)
if m_idx:
m_idx = m_idx if not cpx1_frag1_map else cpx1_frag1_map[m_idx]
lig1_st.atom[m_idx].property[constants.FEP_MAPPING] = atom.index
# get RBs that are common to lig1 and lig2
for rb in analyze.rotatable_bonds_iterator(lig1_st, max_size=MAX_RING_SIZE):
a2, a3 = rb
if not (_valid_mapped_fep_atom(lig1_st.atom[a2]) and
_valid_mapped_fep_atom(lig1_st.atom[a3])):
continue
a2_2 = lig1_st.atom[a2].property[constants.FEP_MAPPING]
a3_2 = lig1_st.atom[a3].property[constants.FEP_MAPPING]
bond_2 = lig2_st.getBond(lig2_st.atom[a2_2], lig2_st.atom[a3_2])
# Check that the rotatable bond is also in lig2_st (DESMOND-10109)
if bond_2 and analyze.is_bond_rotatable(bond_2):
core1_rb.append((a2, a3))
core2_rb.append((a2_2, a3_2))
# get torsions in core for both ligands
core_tors1 = []
core_tors2 = []
for rb1, rb2 in zip(core1_rb[:], core2_rb[:]):
a2_rb1, a3_rb1 = rb1
a1_rb1 = sort_atoms(lig1_st, a2_rb1, a3_rb1, True)
a4_rb1 = sort_atoms(lig1_st, a3_rb1, a2_rb1, True)
if not (a1_rb1 and a4_rb1):
core1_rb.remove(rb1)
core2_rb.remove(rb2)
continue
a2_rb2, a3_rb2 = rb2
a1_rb2 = lig1_st.atom[a1_rb1].property[constants.FEP_MAPPING]
a4_rb2 = lig1_st.atom[a4_rb1].property[constants.FEP_MAPPING]
if not (a1_rb2 and a4_rb2):
core1_rb.remove(rb1)
core2_rb.remove(rb2)
continue
core_tors1.append((a1_rb1, a2_rb1, a3_rb1, a4_rb1))
core_tors2.append((a1_rb2, a2_rb2, a3_rb2, a4_rb2))
lone_tors1 = _getLoneTorsionFEP(lig1_st, core1_rb)
lone_tors2 = _getLoneTorsionFEP(lig2_st, core2_rb)
return core_tors1, core_tors2, lone_tors1, lone_tors2
[docs]def lig_to_fullsys_aid(torsion_list, st):
"""
Convert ligand atom IDs to full-system atom IDs
"""
a, b, c, d = torsion_list
return [
st.atom[a].property[constants.ORIGINAL_INDEX],
st.atom[b].property[constants.ORIGINAL_INDEX],
st.atom[c].property[constants.ORIGINAL_INDEX],
st.atom[d].property[constants.ORIGINAL_INDEX]
]
[docs]def getTorsionKeywords(ligand_st, ligand_asl):
"""
Returns keywords for ligand torsion analysis
:param ligand_st: Structure of the ligand
:type ligand_st: schrodinger.structure.Structure
:param ligand_asl: The ASL to describe the ligand
:type ligand_asl: str
:return: Keyword(s) for calculation
:rtype: sea.Map
"""
new_ligand_st = canonicalize_ligand(ligand_st)
rb_list = get_rotatable_bonds(new_ligand_st)
for rb in rb_list:
ark_kw = sea.Map()
ark_kw["Torsion"] = sea.Map()
ark_kw["Torsion"]["a1"] = rb[0]
ark_kw["Torsion"]["a2"] = rb[1]
ark_kw["Torsion"]["a3"] = rb[2]
ark_kw["Torsion"]["a4"] = rb[3]
ark_kw["Torsion"]["Unit"] = "Degrees"
ark_kw["Torsion"]["Panel"] = 'pl_interact_survey'
ark_kw["Torsion"]["Tab"] = 'l_torsions_tab'
ark_kw["Torsion"]["LigandASL"] = ligand_asl
yield ark_kw
[docs]def getPPIKeywords(protein_asl):
"""
Returns keywords for PLI analysis
:param protein_asl: The ASL to describe the entire protein
:type protein_asl: str
"""
ark_kw = sea.Map()
ark_kw["ProtProtInter"] = sea.Map()
ark_kw["ProtProtInter"]["ASL"] = protein_asl
return ark_kw
[docs]def getPLISKWList(cms_st,
ligand_st,
ligand_asl,
ref_struct_fname,
frame=0,
protein_asl=PROTEIN_ASL,
want_rmsd=True,
want_prmsf=True,
want_lrmsf=True,
want_pli=True,
want_ltorsion=True,
want_lprops=True,
want_ppi=True,
protein_fep=False,
metal_asl=None):
"""
Generate the entire keyword list for all PLI calculations. Also returns
the ligand_asl used in those keywords.
:param ligand_st: Structure of the ligand
:type ligand_st: schrodinger.structure.Structure
:param ligand_asl: The ASL to describe the ligand
:type ligand_asl: str
:param ref_struct_fname: Path to the structure to do the RMSD against
:type ref_struct_fname: None or `str`
:param frame: The frame to take the RMSF/RMSD against, ignored if
ref_struct_fname!=None
:type frame: int
:param protein_asl: The ASL to describe the entire protein
:type protein_asl: str
:param want_rmsd: Whether to add RMSD keywords
:type want_rmsd: bool
:param want_prmsf: Whether to add RMSF keywords for protein
:type want_prmsf: bool
:param want_lrmsf: Whether to add an RMSF ligand keyword
:type want_lrmsf: bool
:param want_pli: Whether to add a PLI keyword
:type want_pli: bool
:param want_ppi: Whether to add a PPI keyword
:type want_ppi: bool
:param want_ltorsion: Whether to add a ligand torsions keyword
:type want_ltorsion: bool
:param want_lprops: Whether to add ligand properties keyword
:type want_lprops: bool
:param protein_fep: Whether the input system is a protein_fep job
:type protein_fep: bool
:param metal_asl: The ASL to describe metals or ions.
If None, use '(ions) or (metals) or (metalloids)'.
:type metal_asl: str or None
:return: Keyword(s) for calculation
:rtype: sea.List
"""
l = sea.List()
lig_asl = '' if protein_fep else ligand_asl
if want_rmsd:
for ark_kw in getRMSDKeywords(protein_asl,
lig_asl,
ref_struct_fname,
frame=frame):
l.append(ark_kw)
if ligand_asl:
for fitby in [None, protein_asl]:
ark_kw = getRMSDLigandKW(ligand_asl,
fitby,
ref_struct_fname,
frame=frame)
l.append(ark_kw)
if want_prmsf:
for ark_kw in getSSEProtKeywords(protein_asl, lig_asl):
l.append(ark_kw)
for ark_kw in getRMSFProtKeywords(protein_asl,
lig_asl,
ref_struct_fname,
frame=frame):
l.append(ark_kw)
if want_lrmsf:
ark_kw = getRMSFLigandKW(ligand_asl + " and not atom.e H",
ligand_asl + " and not atom.e H",
ref_struct_fname)
l.append(ark_kw)
if protein_asl:
ark_kw = getRMSFLigandKW(protein_asl,
ligand_asl + " and not atom.e H",
ref_struct_fname)
l.append(ark_kw)
if want_pli:
l.append(getProtLigInterKW(protein_asl, ligand_asl,
metal_asl=metal_asl))
if want_ppi:
l.append(getPPIKeywords(protein_asl))
if want_ltorsion:
for ark_kw in getTorsionKeywords(ligand_st, ligand_asl):
l.append(ark_kw)
if want_lprops:
# if want_rmsd is not selected, then ligand RMSD will not be
# calculated. So bool 'calcRMSD' is passed for RMSD calculation.
calcRMSD = not want_rmsd
for ark_kw in getLigandPropsKeywords(ligand_asl, calcRMSD,
ref_struct_fname, frame):
l.append(ark_kw)
return l