Source code for schrodinger.application.vsw.vswinput
"""
Class for writing VSW (not general Pipeline) input files.
Copyright Schrodinger, LLC. All rights reserved.
"""
# Contributors: Matvey Adzhigirey, Jon-Eric Eufemio
import os
import sys
import schrodinger.job.util as jobutil # To check for Epik installation
import schrodinger.pipeline.input as vswinput
from schrodinger.structutils import structalign
# Check whether SCHRODINGER_PYTHON_DEBUG is set for debugging:
from schrodinger.utils import cmdline
from schrodinger.utils import log
from schrodinger.utils.ligfilter import FILTERFILE_EXT
DEBUG = (log.get_environ_log_level() <= log.DEBUG)
logger = log.get_output_logger("vswinput")
VSW_COMPOUNDFIELD = 's_vsw_compound_code'
GPHDB_COMPOUNDFIELD = 's_m_title' # Ev:105243
VARIANTFIELD = 's_vsw_variant'
PHASEDB = "PhaseDB"
STRUCTURES = "Structures"
GRID = "Grid"
[docs]def parse_command_line():
"""
Returns the path to input file, as determined from command line option parser.
"""
usage = """
To specify a parameters file that will be read to set the job parameters in the VSW job input file:
$SCHRODINGER/run vsw_input.py -p <job parameter file>"""
description = """Description: This script generates a VSW input file. The parameters used for the VSW input file can be specified in a parameters file."""
parser = cmdline.SingleDashOptionParser(usage=usage,
description=description)
parser.add_option("-p",
"-parameters",
type="string",
dest="job_params",
default="",
metavar="FILE",
help="Job parameters file")
(options, args) = parser.parse_args()
return options.job_params
[docs]class VSWWriter():
"""
Use writeInputFile() method
"""
[docs] def __init__(self, jp=None):
"""
Initialize from parameters dictionary.
"""
# Initialize to defaults:
# minimization method
self.prepare = False
# Controls for the MergeDuplicates stage:
self.merge_duplicates = False
self.desalt = False
self.neutralize = False
self.merge_props = False
self.dockligs_files = None
self.dockligs_db = None
self.dockligs_dbsubset = None
# Triggers phasedb_manage -unique smiles option.
self.remove_duplicates = False
# If True, the <db>_unique file will be deleted before adding ligands.
self.delete_unique_file = False
self.regularize = False
self.ph = 7.0
self.pht = 2.0
self.useepik = True
self.metal_binding = False
self.mixligs = True
self.nringconfs = 1
self.stereo_source = 'geometry'
self.uniquefield = "NONE"
# NONE means to retitle by position
self.recombine = True
self.hiprot = False
self.sample_rings = False
self.maxstereo = 4
self.lipinski = False
self.reactive = False
self.runqikprop = False
self.precustomfilter = None # Will only run IF preparing
self.postcustomfilter = None # Will run before docking
self.titlefilterfile = None
self.receptors = []
self.align_receptors = False
self.lvdw = 0.8
self.lig_ccut = 0.15
self.maxatoms = 300
self.maxrotbonds = 50
self.res_inter_radius = None # <None> means disable by default
self.runhtvs = False
# whether to keep number of structures as opposed to percentage
self.htvs_donum = False
self.htvsper = 10.0
self.htvsnum = 1
self.htvsposes = 1
self.htvskeep = 'All states'
self.htvsconfs = 'confgen'
self.htvsusecons = False
self.htvsamide = 'penal'
self.htvs_postdock = True
self.htvs_postdockstrain = False
self.htvs_forcefield = None
self.htvs_epik_penalties = False
self.htvs_enhanced_planarity = False # Ev:120767
self.runsp = False
self.sp_donum = False
self.spper = 10.0
self.spnum = 1
self.spposes = 1
self.spkeep = 'All states'
self.spconfs = 'confgen'
self.spusecons = False
self.spamide = 'penal'
self.spwritedesc = False
self.sp_nsampling = 1
self.sp_postdock = True
self.sp_postdockstrain = False
self.sp_forcefield = None
self.sp_epik_penalties = False
self.sp_enhanced_planarity = False # Ev:120767
self.runxp = False
self.xpgenconfs = False
self.xpwritedesc = False
self.xpcompmax = False
self.xpcompmax_cons = True
self.mmgbsa = False
self.xp_donum = False
self.xpper = 100.0
self.xpnum = 1
self.xpposes = 1
self.xpkeep = 'All states'
self.xpconfs = 'confgen'
self.xpusecons = False
self.xpamide = 'penal'
self.xp_postdock = True
self.xp_postdockstrain = False
self.xp_forcefield = None
self.xp_epik_penalties = False
self.xp_enhanced_planarity = False # Ev:120767
self.phasedb_path = None
self.phasedb_new = True
self.phasedb_confs_mode = "auto"
self.phasedb_maxconfs = 100
self.phasedb_maxperbond = 10
self.phasedb_amide = "vary"
self.phasedb_energy_window = 10.0
self.phasedb_sample_method = "rapid"
self.jobname = 'vsw'
# Where to write the input file:
self.input_file_path = None # None: <jobname>.inp
# Adjust to specified parameters dict:
if jp:
self.__dict__.update(jp)
[docs] def writeDockingStage(self,
mode,
receptor,
input,
output,
outtype,
gencodes_out_set=None,
frags=False,
recombine=True,
unique_field=None,
multiple_confs=False):
"""
Writes stage that describes how the docking will be executed.
"""
if receptor.identifier == '':
receptor_suffix = ''
else:
receptor_suffix = '_' + receptor.identifier
if frags:
sname = "DOCK_%s%s_FRAGS" % (mode, receptor_suffix)
else:
sname = "DOCK_%s%s" % (mode, receptor_suffix)
sclass = "glide.DockingStage"
ingrid = "GRID%s" % receptor_suffix
#outputs = [output]
lmode = mode.lower()
precision = mode.upper()
docking_uniquefield = None
if recombine:
# Determine good min/max num of ligand per subjobs:
if frags:
method = "confgen"
else:
method = getattr(self, '%sconfs' % lmode)
# Determine the approximate size of a 1 hour job:
min_job_size = None
if precision == "HTVS":
if method == "inplace" or method == "refineinput":
msg = "ERROR: %s unavailable for %s" % (method, precision)
raise RuntimeError(msg)
elif method == "rigid":
min_job_size = 8000
elif method == "confgen":
min_job_size = 4000
elif precision == "SP":
if method == "inplace":
min_job_size = 10000 # 25,000 ~ 1 hour; 10,000 ~ 25min
elif method == "refineinput":
min_job_size = 5000 # 8,000 ~ 1 hour; 5,000 ~ 40min
elif method == "rigid":
min_job_size = 2000
elif method == "confgen":
min_job_size = 300
elif precision == "XP":
if method == "inplace":
min_job_size = 4000
elif method == "refineinput":
min_job_size = 100
elif method == "rigid":
min_job_size = 200
elif method == "confgen":
min_job_size = 20
else:
# Should never happen
msg = 'ERROR: Invalid precision mode: "%s"' % precision
raise RuntimeError(msg)
if min_job_size is None:
# Should no longer ever happen
msg = 'ERROR: docking method "%s" is unavailable for precision "%s"' % (
method, precision)
raise RuntimeError(msg)
max_job_size = min_job_size * 10 # about 10 hours
kw = {
# Use -NJOBS to determine number of subjobs
'NUMOUT': 'njobs',
'OUTFORMAT': 'maegz',
'MIN_SUBJOB_STS': min_job_size,
'MAX_SUBJOB_STS': max_job_size,
}
recom_output = sname + '_INPUT'
if gencodes_out_set:
# If generating variant codes:
kw.update({
'GENCODES': True,
'OUTCOMPOUNDFIELD': VSW_COMPOUNDFIELD,
'OUTVARIANTFIELD': VARIANTFIELD,
# This is the first PRE_DOCK stage, generate codes from the
# user-specified uniquefield:
'UNIQUEFIELD': unique_field,
})
docking_uniquefield = VSW_COMPOUNDFIELD
# Only 1 stage should have RECOMBINE_OUT as output
recom_output = gencodes_out_set
elif frags:
# Docking fragments (they don't have variant codes
kw.update({
'GENCODES': True,
'OUTCOMPOUNDFIELD': VSW_COMPOUNDFIELD,
'OUTVARIANTFIELD': VARIANTFIELD,
'UNIQUEFIELD': "s_m_title",
})
docking_uniquefield = "s_m_title"
else:
# Not generating variant codes:
kw['GENCODES'] = False # Ev:124913
if frags:
# Docking fragments, they don't have the user-specified
# unique property:
kw['UNIQUEFIELD'] = 's_m_title'
docking_uniquefield = "s_m_title"
else:
kw['UNIQUEFIELD'] = unique_field
docking_uniquefield = unique_field
recomsname = 'PRE_' + sname
recomsclass = 'gencodes.RecombineStage'
self.writer.addStage(recomsname, recomsclass, [input],
[recom_output], kw)
input = recom_output
else:
docking_uniquefield = unique_field
# Should never happen:
if docking_uniquefield is None:
raise RuntimeError("unique-field logic error")
# Beginning of the docking stage keywords
kw = {
# We are now always including RecombineStage if recombining is
# necessary
'RECOMBINE': False,
'PRECISION': mode,
'UNIQUEFIELD': docking_uniquefield,
}
# If docking fragments for XPVis (Ev:63790) keep all of them:
if frags:
kw['PERCENT_TO_KEEP'] = 100
kw['DOCKING_METHOD'] = 'confgen'
else:
if getattr(self, '%s_donum' % lmode):
kw['NUM_TO_KEEP'] = getattr(self, '%snum' % lmode)
else:
kw['PERCENT_TO_KEEP'] = getattr(self, '%sper' % lmode)
kw['DOCKING_METHOD'] = getattr(self, '%sconfs' % lmode)
if frags: # For fragments, keep 9 poses per ligand (Ev:68950)
poses_per_lig = 9
else:
poses_per_lig = getattr(self, '%sposes' % lmode)
kw['POSES_PER_LIG'] = poses_per_lig
if mode == 'SP':
kw['WRITE_XP_DESC'] = self.spwritedesc
kw['NENHANCED_SAMPLING'] = self.sp_nsampling
elif mode == 'XP':
kw['WRITE_XP_DESC'] = self.xpwritedesc
keepmatch = {
'Only best scoring state': 'YES',
'Only best scoring pose': 'YES',
'All good scoring states': 'NO',
'All good scoring poses': 'NO',
'All states': 'NO',
'best': 'YES', # NEW
'good': 'NO', # NEW
'all': 'NO', # NEW
}
keep_best = keepmatch[getattr(self, '%skeep' % lmode)]
kw.update({
'BEST_BY_TITLE': keep_best,
'LIG_VSCALE': self.lvdw,
'LIG_CCUT': self.lig_ccut,
'MAXATOMS': self.maxatoms,
'MAXROTBONDS': self.maxrotbonds,
})
# Whether to report per-residue interactions VSW-824:
res_inter_radius = self.res_inter_radius
if res_inter_radius:
kw['WRITE_RES_INTERACTION'] = True
kw['RADIUS_RES_INTERACTION'] = res_inter_radius
# Match pull down output with the short name to pass to MMIM:
amide_mode = getattr(self, '%samide' % lmode)
# If any conf other than confgen (flexible) is selected, is it
# necessary to set the keyword to default value?????
# NO for SP and YES for XP
forcefield = getattr(self, '%s_forcefield' % lmode)
postdock = getattr(self, '%s_postdock' % lmode)
postdockstrain = getattr(self, '%s_postdockstrain' % lmode)
kw['AMIDE_MODE'] = amide_mode
kw['POSE_OUTTYPE'] = outtype
if forcefield is not None:
kw['FORCEFIELD'] = forcefield
kw.update({
'POSTDOCK': postdock,
'POSTDOCKSTRAIN': postdockstrain, # Ev:70957
'COMPRESS_POSES': True, # Ev:87161
})
# Ev:80478 & Ev:90171
epik_penalties = getattr(self, '%s_epik_penalties' % lmode)
kw['EPIK_PENALTIES'] = epik_penalties
enhanced_planarity = getattr(self, '%s_enhanced_planarity' % lmode)
kw['FORCEPLANAR'] = enhanced_planarity
if multiple_confs:
# If input ligands have multiple conformers, we need to make
# sure we do not canonicalize them into a single conformer:
kw['CANONICALIZE'] = False
self.htvs_enhanced_planarity = False # Ev:120767
if not frags or (frags and self.xpcompmax_cons):
# Ev:72134 Setup constraints for fragments only if requested
# Setup HBOND and METAL constraints:
if getattr(self, '%susecons' % lmode): # Use constraints:
# At least one specified for this receptor
if receptor.usecons:
if receptor.features:
# If at least one positional constraint is present
# (requires features):
feature_dict = {}
for feature_num, feature in receptor.features.items():
feature_config = {}
for pattern_num, pattern in enumerate(
feature.patterns, start=1):
pattern_name = "PATTERN%i" % pattern_num
atoms_str = ",".join(map(str, pattern.atoms))
if pattern.exclude:
inc_exc_str = "exclude"
else:
inc_exc_str = "include"
feature_config[pattern_name] = "%s %s %s" % (
pattern.smarts, atoms_str, inc_exc_str)
kw["FEATURE:%i" % feature_num] = feature_config
# Write a single constraint group to the VSW input file:
kw['CONSTRAINT_GROUP:1'] = {
"USE_CONS": receptor.usecons,
"NREQUIRED_CONS": receptor.nrequired,
}
# Setup CORE constraint:
if receptor.usecore: # Core defined for this receptor
core_file = '%s-receptor%s-core.mae' % \
(self.jobname, receptor_suffix)
receptor.core_structure.write(core_file)
kw.update({
'CORE_RESTRAIN': True,
'USE_REF_LIGAND': True,
'REF_LIGAND_FILE': core_file,
'CORE_POS_MAX_RMSD': receptor.core_tolerance,
'CORE_DEFINITION': receptor.core_atoms,
'CORE_FILTER': receptor.core_skipnonmatching,
})
if receptor.core_atoms == 'smarts':
# Which atoms of the ligand define the core:
kw['CORE_SMARTS'] = receptor.core_smarts
kw['CORE_ATOMS'] = receptor.core_atom_list
self.writer.addStage(sname, sclass, [input, ingrid], [output], kw)
[docs] def writeGenConfsStage(self, receptor_suffix, input, output):
"""
Writes stages that generate MMFFs and OPLS conformers from the input,
combine them with the input, and put into output.
"""
sname = "MMFFS" + receptor_suffix
sclass = "macromodel.ConfSearchStage"
out1 = "MMFFS_OUT" + receptor_suffix
kw = {
'FORCE_FIELD': "MMFFs",
'MINI_METHOD': "TNCG",
'SOLVENT': "Water",
'CUTOFF': "Extended",
}
self.writer.addStage(sname, sclass, [input], [out1], kw)
sname = "OPLS" + receptor_suffix
sclass = "macromodel.ConfSearchStage"
out2 = "OPLS_OUT" + receptor_suffix
kw = {
'FORCE_FIELD': "OPLS_2005",
'MINI_METHOD': "TNCG",
'SOLVENT': "Water",
'CUTOFF': "Extended",
}
self.writer.addStage(sname, sclass, [input], [out2], kw)
sname = "COMBINE" + receptor_suffix
sclass = "combine.CombineStage"
kw = {
'LABELFIELD': "s_vsw_conformer_field",
'LABELS': ["Original", "MMFFS", "OPLS"],
}
self.writer.addStage(sname, sclass, [input, out1, out2], [output], kw)
[docs] def writePreFilterStage(self, input, output, custom_filter_text):
"""
Writes the stage that filters by specified LigFilter criteria
"""
custom_filter_file = '%s-custom-filter.%s' % \
(self.jobname, FILTERFILE_EXT)
fh = open(custom_filter_file, 'w')
fh.write(custom_filter_text)
fh.close()
sname = "CUSTOMFILTER"
sclass = "filtering.LigFilterStage"
kw = {'FILTER_FILE': custom_filter_file}
self.writer.addStage(sname, sclass, [input], [output], kw)
[docs] def writePrepareStages(self, lastoutput, recombine):
if self.precustomfilter:
output = 'CUSTOMFILTER_OUT'
self.writePreFilterStage(lastoutput, output, self.precustomfilter)
lastoutput = output
sname = 'LIGPREP'
sclass = "ligprep.LigPrepStage"
input = lastoutput
lastoutput = 'LIGPREP_OUT'
# This stage is going to be the first in the workflow, so set the compound code,
# unless the user chose to not recombine:
if recombine:
compoundfield = VSW_COMPOUNDFIELD
kw = {
'RECOMBINE': True,
'RETITLE': True,
'MIXLIGS': self.mixligs,
'SKIP_BAD_LIGANDS': True,
'UNIQUEFIELD': self.uniquefield,
'OUTCOMPOUNDFIELD': compoundfield,
}
# For the docking stages:
self.uniquefield = VSW_COMPOUNDFIELD
else:
kw = {
'RECOMBINE': False,
'RETITLE': False,
'UNIQUEFIELD': self.uniquefield,
}
kw.update({
'USE_EPIK': self.useepik,
'METAL_BINDING': self.metal_binding,
'PH': self.ph,
'PHT': self.pht,
'NRINGCONFS': self.nringconfs,
'COMBINEOUTS': False,
'STEREO_SOURCE': self.stereo_source,
'NUM_STEREOISOMERS': 32, # For -s option
'MAX_STEREOISOMERS': self.maxstereo, # For -m option
})
if self.merge_duplicates:
# If MergeDuplicatesStage is run, the structures are already
# regularized:
self.regularize = False
kw['REGULARIZE'] = self.regularize
self.writer.addStage(sname, sclass, [input], [lastoutput], kw)
if self.sample_rings:
output = 'SAMPLERINGS_OUT'
keywords = {'OUTCONFS_PER_SEARCH': 10} # Ev:91804
self.writer.addStage('SAMPLERINGS', 'macromodel.SampleRingsStage',
[lastoutput], [output], keywords)
lastoutput = output
# PostLigPrep stage:
if recombine:
uniquefield = VSW_COMPOUNDFIELD
else:
if self.uniquefield == "NONE":
uniquefield = "s_m_title" # Generated by LigPrep
else:
uniquefield = self.uniquefield
keywords = {
'UNIQUEFIELD': uniquefield,
'OUTVARIANTFIELD': VARIANTFIELD,
# creates same number of QikProp jobs as LigPrep
'PRESERVE_NJOBS': True,
}
keywords['REMOVE_PENALIZED_STATES'] = 'YES' if self.hiprot else 'NO'
output = 'POSTLIGPREP_OUT'
self.writer.addStage('POSTLIGPREP', 'ligprep.PostLigPrepStage',
[lastoutput], [output], keywords)
lastoutput = output
return lastoutput
[docs] def writeInputFile(self):
"""
Writes the VSW input file
Raises a RuntimeError on error.
"""
if self.recombine:
compoundfield = VSW_COMPOUNDFIELD
else:
compoundfield = self.uniquefield
pullmatch = {
'Only best scoring state': VARIANTFIELD,
'Only best scoring pose': VARIANTFIELD,
'All good scoring states': VARIANTFIELD,
'All good scoring poses': VARIANTFIELD,
'All states': compoundfield,
'best': VARIANTFIELD, # NEW
'good': VARIANTFIELD, # NEW
'all': compoundfield, # NEW
}
# Used by the PULL stage:
htvspull_prop = pullmatch[self.htvskeep]
sppull_prop = pullmatch[self.spkeep]
if self.input_file_path:
vsw_input_file = self.input_file_path
else:
vsw_input_file = self.jobname + '.inp'
# Ev:125529
txt = "########## Virtual Screening Workflow Input File ###############\n\n"
txt += "# Run as: $SCHRODINGER/vsw <inputfile> \n"
if False: # JON_TEST
txt += "\n"
txt += "\n"
txt += "# Print sorted self.__dict__ key/value pairs): \n"
txt += "# \n"
jp_keys = list(self.__dict__)
sorted_jp_keys = sorted(jp_keys)
jp_index = 1
for key in sorted_jp_keys:
txt += "# " + str(jp_index).rjust(2) + " " + str(key).rjust(20) + " " + \
str(type(self.__dict__[key])).ljust(
25) + " " + str(self.__dict__[key]) + "\n"
jp_index += 1
txt += "\n"
txt += "\n"
txt += "# BEGIN Print jp as assignments: \n"
jp_index = 1
for key in sorted_jp_keys:
jp_keyname = "self." + str(key)
variable_quoting = ''
if type(self.__dict__[key]).__name__ == 'str':
variable_quoting = '"'
jp_value = variable_quoting + \
str(self.__dict__[key]) + variable_quoting
txt += "# " + jp_keyname.ljust(20) + " = " + jp_value.ljust(50) + \
" # " + \
str(type(self.__dict__[key])) + \
"\n"
txt += "# END Print jp as assignments: \n"
txt += "\n"
self.writer = vswinput.Writer(vsw_input_file, comment=txt)
logger.debug("ADDING VARIABLES...")
if self.dockligs_files is not None:
# Set input ligands:
self.writer.addVar('ORIGINAL_LIGANDS', STRUCTURES,
self.dockligs_files)
lastoutput = "ORIGINAL_LIGANDS"
elif self.dockligs_db is not None:
# Set input database:
dbvar = 'INPUT_DATABASE'
self.writer.addVar(dbvar, PHASEDB, self.dockligs_db)
lastoutput = 'EXPORTED_LIGANDS'
keywords = {'OUTFORMAT': 'maegz', 'MAX_CONFS': 1}
# Save the full path to the file (Ev:127191)
subsetfile = self.dockligs_dbsubset
if subsetfile:
keywords['SUBSET_FILE'] = subsetfile
self.writer.addStage('EXPORT', 'phase.DBExportStage', [dbvar],
[lastoutput], keywords)
else:
raise ValueError("Input ligands are not specified")
# Make sure that no docking takes place if unchecked:
if not len(self.receptors) > 0:
self.runhtvs = False
self.runsp = False
self.runxp = False
if self.align_receptors:
# Ev:81292 Align all receptors (except pre-generated grids) to the
# first receptor:
reference_st = None
mobile_sts = []
for rec in self.receptors:
if rec.generate_grid:
if not reference_st:
reference_st = rec.model
else:
mobile_sts.append(rec.model)
if not reference_st:
raise RuntimeError(
"Can't align receptors because no receptors are specified that are not pre-generated grids."
)
if not mobile_sts:
raise RuntimeError(
"Can't align receptors because only one receptor is specified (excluding pre-generated grids)."
)
# Will align the actual CTs stored in the Receptor objects:
ska = structalign.StructAlign()
ska.align(reference_st, mobile_sts)
# Set input grids:
for rec in self.receptors:
if rec.identifier == '':
receptor_suffix = ''
else:
receptor_suffix = '_' + rec.identifier
if not rec.generate_grid:
gridvarname = 'GRID' + receptor_suffix
self.writer.addVar(gridvarname, GRID, rec.path)
else:
# Write the receptor (used to be in Workspace) to a file:
model_file = '%s-receptor%s.mae' % \
(self.jobname, receptor_suffix)
rec.model.write(model_file)
modelvarname = 'MODEL' + receptor_suffix
self.writer.addVar(modelvarname, STRUCTURES, [model_file])
# Set fragments library (if calculating maximum values Ev:63790):
fragsvar = None
if self.runxp and self.xpcompmax:
glide_exec = jobutil.hunt('glide')
if not glide_exec:
msg = "ERROR: Glide is not installed. You can still submit the job to a remote host that has Glide installed; but only if you disable the option to compute maximum values by docking fragments."
raise RuntimeError(msg)
glide_dir = os.path.split(os.path.split(
jobutil.hunt('glide'))[0])[0]
fragments_file = os.path.join(glide_dir, 'data', 'glide',
'opt-xp-50frags_epik.mae.gz')
if not os.path.isfile(fragments_file):
msg = "Could not locate the fragments file: %s" % fragments_file
raise RuntimeError(msg)
fragsvar = 'FRAGMENTS'
self.writer.addVar(fragsvar, STRUCTURES, [fragments_file])
# Never convert to SMILES if not preparing ligands:
if not self.prepare and self.regularize:
print(
'WARNING: vswinput: prepare is %s and regularize is %s. Will not regularize'
% (self.prepare, self.regularize))
self.regularize = False
if self.merge_duplicates:
output = 'NODUPLICATES'
kw = {
'NEUTRALIZE': self.neutralize,
'DESALT': self.desalt,
'MERGE_PROPS': self.merge_props,
'OUTFORMAT': 'sdf',
'SMILES_FIELD': 'VendorSMILES',
'CODE_FIELD': None,
'CODE_PREFIX': None,
'REQUIRE_PROPS': False,
}
self.writer.addStage('MERGE_DUPLICATES',
'filtering.MergeDuplicatesStage', [lastoutput],
[output], kw)
lastoutput = output
if self.titlefilterfile:
# Write stage for filtering based on title (or other prop):
# (Supports SMILES files as well)
output = 'TITLEFILTER_OUT'
kw = {
'PROPERTY': self.uniquefield,
'VALUE_FILE': self.titlefilterfile,
}
self.writer.addStage('TITLEFILTER', 'filtering.SubSetStage',
[lastoutput], [output], kw)
lastoutput = output
if self.prepare:
logger.debug("ADDING PREPARE STAGES")
lastoutput = self.writePrepareStages(lastoutput, self.recombine)
else: # not self.prepare
pass
logger.debug("ADDING PRE-DOCKING STAGES")
# Gets here whether or not Prepare ligands is checked
if self.runqikprop:
recombine_ligands = (not self.prepare and self.recombine)
kw = {'RECOMBINE': recombine_ligands}
output = "QIKPROP_OUT"
self.writer.addStage('QIKPROP', 'qikprop.QikPropStage',
[lastoutput], [output], kw)
lastoutput = output
if self.lipinski:
sname = 'LIPINSKI_FILTER'
sclass = 'filtering.LigFilterStage'
input = lastoutput
lastoutput = 'LIPFILTER_OUT'
conditions = [
'r_qp_mol_MW <= 500', 'r_qp_QPlogPo/w <= 5',
'r_qp_donorHB <= 5', 'r_qp_accptHB <= 10'
]
kw = {'CONDITIONS': conditions}
self.writer.addStage(sname, sclass, [input], [lastoutput], kw)
if self.reactive:
sname = 'REACTIVE_FILTER'
sclass = 'filtering.LigFilterStage'
input = lastoutput
lastoutput = 'REACTIVEFILTER_OUT'
conditions = ['Reactive_groups == 0']
kw = {'CONDITIONS': conditions}
self.writer.addStage(sname, sclass, [input], [lastoutput], kw)
if self.postcustomfilter:
output = 'CUSTOMFILTER_OUT'
self.writePreFilterStage(lastoutput, output, self.postcustomfilter)
lastoutput = output
if self.xpgenconfs and self.runxp and not self.runhtvs and not self.runsp:
# Ev:53413 Generate conformers once for all grids if only XP
# docking is used:
output = "MIN_COMBINED"
self.writeGenConfsStage("", lastoutput, output)
lastoutput = output
self.prepared_ligand_set = lastoutput
if self.prepare or not self.receptors:
# Ev:131515 Return the last ligand set to the user if
# preparation took place OR no docking took place
self.writer.userOutput(self.prepared_ligand_set)
if not self.receptors:
# We are not performing docking
self.writer.setStructureOutput(self.prepared_ligand_set)
# Do not use VARIANTFIELD for pulling if codes were not generated.
# This would happen if recombine is False and not Preparing.
if (not self.prepare) and (not self.recombine):
# FIXME: Would be better to check if the input files have
# the variant property in them...
if self.uniquefield == 'NONE':
htvspull_prop = 's_m_title'
sppull_prop = 's_m_title'
else:
htvspull_prop = self.uniquefield
sppull_prop = self.uniquefield
docking_outputs = []
docking_offsets = []
docking_identifiers = []
# Add docking stages for all receptors:
logger.debug("ADDING DOCKING STAGES")
for rec in self.receptors:
if rec.identifier == '':
receptor_suffix = ''
else:
receptor_suffix = '_' + rec.identifier
# If recombining (and ligands were not prepared), use a special
# pull set for this receptor:
if not self.prepare and self.recombine:
pull_from_set = 'RECOMBINE_OUT' + receptor_suffix
# First docking stage will produce this:
gencodes_out_set = pull_from_set
else:
# Use prepared ligands (or original ligands):
# (or not preparing and not recombining)
pull_from_set = self.prepared_ligand_set
gencodes_out_set = None
gridvarname = 'GRID' + receptor_suffix
if rec.generate_grid:
# Generate the grids:
modelvarname = 'MODEL' + receptor_suffix
sname = "GRIDGEN" + receptor_suffix
sclass = "glide.GridgenStage"
inner_box_length = 10 # center of ligand box
kw = {
'GRID_CENTER': rec.box_center, # Ev:100860
# box length (advanced options; int)
'INNERBOX': inner_box_length,
# box length + ligsize
'OUTERBOX': float(inner_box_length) + rec.lig_size,
'RECEP_VSCALE': rec.rvdw,
'WRITEZIP': True,
}
# Add constraints to be generated:
hc = ['%s' % c for c in rec.hbond_constraints]
if hc:
# hc is a list of "label atomnum" strings
kw['HBOND_CONSTRAINTS'] = hc
mc = ['%s' % c for c in rec.metal_constraints]
if mc:
# mc is a list of "label atomnum" strings
kw['METAL_CONSTRAINTS'] = mc
pc = ['%s' % c for c in rec.posit_constraints]
if pc:
# pc is a list of "label x y z radius" strings
kw['POSIT_CONSTRAINTS'] = pc
nc = ['%s' % c for c in rec.noe_constraints]
if nc:
# pc is a list of "label x y z rmin rmax" strings
kw['NOE_CONSTRAINTS'] = nc
self.writer.addStage(sname, sclass, [modelvarname],
[gridvarname], kw)
lastoutput = self.prepared_ligand_set
# HTVS output after being pulled from the input set (optimization):
htvs_out_orig = None
if self.runhtvs:
first_docking = True
last_docking = (not self.runsp and not self.runxp)
# Will NOT recombine ONLY when this is the first stage AND
# self.recombine is False:
recombine = self.prepare or self.recombine or (
not first_docking)
if not last_docking:
outtype = 'LIB'
else:
outtype = 'PV'
# This is the first docking stage. Honor the users's
# uniquefield:
stage_unique_field = self.uniquefield
output = "HTVS_OUT" + receptor_suffix
self.writeDockingStage('HTVS',
rec,
lastoutput,
output,
outtype,
gencodes_out_set,
recombine=recombine,
unique_field=stage_unique_field)
gencodes_out_set = None
lastoutput = output
pull = False
if self.runsp: # Next stage is SP
if self.spconfs not in ['inplace', 'refineinput']:
pull = True
elif self.runxp: # Next stage is XP
if self.xpconfs not in ['inplace', 'refineinput']:
pull = True
if pull:
sname = "PULL_HTVS" + receptor_suffix
sclass = "pull.PullStage"
in1 = lastoutput
in2 = pull_from_set
lastoutput = "HTVS_OUT_ORIG" + receptor_suffix
kw = {'UNIQUEFIELD': htvspull_prop}
self.writer.addStage(sname, sclass, [in1, in2],
[lastoutput], kw)
htvs_out_orig = lastoutput
if self.runsp:
first_docking = (not self.runhtvs)
last_docking = (not self.runxp)
# Will NOT recombine ONLY when this is the first stage AND
# self.recombine is False:
recombine = self.prepare or self.recombine or (
not first_docking)
if not last_docking:
outtype = 'LIB'
else:
outtype = 'PV'
# If this is the first docking stage, honor the user's uniquefield.
# If the user chose to not recombine, also honor the user's
# uniquefield.
if first_docking or not self.recombine:
stage_unique_field = self.uniquefield
else:
stage_unique_field = VSW_COMPOUNDFIELD
output = "SP_OUT" + receptor_suffix
self.writeDockingStage('SP',
rec,
lastoutput,
output,
outtype,
gencodes_out_set,
recombine=recombine,
unique_field=stage_unique_field)
gencodes_out_set = None
lastoutput = output
if self.runxp and self.xpconfs not in [
'inplace', 'refineinput'
]:
sname = "PULL_SP" + receptor_suffix
sclass = "pull.PullStage"
in1 = "SP_OUT" + receptor_suffix
if htvs_out_orig is not None:
# HTVS output pulled from the original set
# (optimization):
in2 = htvs_out_orig
else:
in2 = pull_from_set
lastoutput = "SP_OUT_ORIG" + receptor_suffix
kw = {'UNIQUEFIELD': sppull_prop}
self.writer.addStage(sname, sclass, [in1, in2],
[lastoutput], kw)
if self.runxp:
first_docking = (not self.runhtvs and not self.runsp)
last_docking = True
# Will NOT recombine ONLY when this is the first stage AND
# self.recombine is False:
recombine = self.prepare or self.recombine or (
not first_docking)
if self.xpgenconfs and (self.runhtvs or self.runsp):
# Ev:53413 Generate conformers before the XP stage if there
# are docking stages before this one:
output = "MIN_COMBINED" + receptor_suffix
self.writeGenConfsStage(receptor_suffix, lastoutput, output)
lastoutput = output
# If this is the first docking stage, honor the user's uniquefield.
# If the user chose to not recombine, also honor the user's
# uniquefield.
if first_docking or not self.recombine:
stage_unique_field = self.uniquefield
else:
stage_unique_field = VSW_COMPOUNDFIELD
output = "XP_OUT%s" % receptor_suffix
self.writeDockingStage('XP',
rec,
lastoutput,
output,
'PV',
gencodes_out_set,
recombine=recombine,
unique_field=stage_unique_field,
multiple_confs=self.xpgenconfs)
gencodes_out_set = None
self.writer.userOutput(output) # Ev:72832
lastoutput = output
if self.xpcompmax:
output = 'XP_OUT%s_xpfrags' % receptor_suffix
# fragsvar was set earlier
stage_unique_field = "s_m_title"
self.writeDockingStage('XP',
rec,
fragsvar,
output,
'PV',
frags=True,
recombine=recombine,
unique_field=stage_unique_field)
self.writer.userOutput(output)
if self.mmgbsa:
output = 'MMGBSA' + receptor_suffix
self.writer.addStage('MMGBSA' + receptor_suffix,
'prime.MMGBSAStage', [lastoutput],
[output])
lastoutput = output
# Output from LAST docking stage
docking_outputs.append(lastoutput)
self.writer.userOutput(lastoutput)
docking_offsets.append(str(rec.offset)) # Offset for this receptor
# Marker for this receptor
docking_identifiers.append(rec.identifier)
logger.debug("ADDING POST-DOCKING STAGES")
if len(docking_outputs) > 1: # Ensemble docking was done
# Merge the docking results using the Merge stage:
# Also keeps 1 (or max_per_lig) conformers per compound
# so needs to run if XP has multiple input structures per root.
final_dock_out = "FINAL_DOCK_OUT"
# If not recombining, the output of the first docking stage won't
# have the VSW_COMPOUNDFIELD property:
if not self.recombine:
# "Use input files directly" radio selected
stage_unique_field = self.uniquefield
else:
stage_unique_field = VSW_COMPOUNDFIELD
kw = {
'MAXPERLIG': 1,
'NREPORT': 10000,
'UNIQUEFIELD': stage_unique_field,
'OFFSETS': docking_offsets,
'MARKERS': docking_identifiers,
}
self.writer.addStage("DOCKMERGE", "glide.MergeStage",
docking_outputs, [final_dock_out], kw)
# Add merged output to useroutput:
self.writer.userOutput(final_dock_out)
self.writer.setStructureOutput(final_dock_out)
elif len(docking_outputs) == 1:
self.writer.setStructureOutput(docking_outputs[0])
logger.debug('WRITING THE INPUT FILE')
self.writer.write()
logger.debug('DONE')
return vsw_input_file
if __name__ == '__main__':
params_file = parse_command_line()
if not params_file:
sys.exit(1)
# check that job parameters file was specified; if it was
# check that it exists before trying to use it.
if not os.path.isfile(params_file):
print('\n Error: Job parameters file "%s" does not exist.' %
params_file)
sys.exit(1)
class vsw_receptor:
pass
receptor = vsw_receptor()
receptor.generate_grid = False
jp = {'receptors': [receptor]}
for param_file_line in open(params_file):
exec(param_file_line)
inputfile = VSWWriter(jp).writeInputFile()
print('Written input file:', inputfile)
# EOF