Skip to content

Commit

Permalink
added creating box around ligand
Browse files Browse the repository at this point in the history
  • Loading branch information
@mattholc
mattholc committed Jul 13, 2023
1 parent b8b4e05 commit efb14bc
Show file tree
Hide file tree
Showing 3 changed files with 91 additions and 47 deletions.
38 changes: 34 additions & 4 deletions meeko/gridbox.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ def box_to_pdb_string(box_center, npts, spacing=0.375):
|5....|./6
|.____|/
1 2
"""

step_x = int(npts[0] / 2.0) * spacing
Expand All @@ -77,7 +77,7 @@ def box_to_pdb_string(box_center, npts, spacing=0.375):
corners.append([center_x + step_x, center_y - step_y, center_z + step_z] ) # 6
corners.append([center_x + step_x, center_y + step_y, center_z + step_z] ) # 7
corners.append([center_x - step_x, center_y + step_y, center_z + step_z] ) # 8

count = 1
res = "BOX"
chain = "X"
Expand All @@ -89,10 +89,10 @@ def box_to_pdb_string(box_center, npts, spacing=0.375):
z = corners[idx][2]
pdb_out += line % (count, "Ne", res, chain, idx+1, x, y, z, "Ne")
count += 1

# center
pdb_out += line % (count+1, "Xe", res, chain, idx+1, center_x, center_y, center_z, "Xe")

pdb_out += "CONECT 1 2" + os_linesep
pdb_out += "CONECT 1 4" + os_linesep
pdb_out += "CONECT 1 5" + os_linesep
Expand Down Expand Up @@ -122,5 +122,35 @@ def is_point_outside_box(point, center, npts, spacing=0.375):
is_outside |= z >= maxcorner[2] or z <= mincorner[2]
return is_outside

def calc_box(fname, padding):
"""Crude PDBQT parsing is used here because it allows input of autosite pdbqts"""
padding = float(padding)
x_min = float('inf')
x_min = float('inf')
y_min = float('inf')
z_min = float('inf')
x_max = float('-inf')
y_max = float('-inf')
z_max = float('-inf')
with open(fname) as f:
for line in f:
if line.startswith('ATOM') or line.startswith('HETATM'):
x = float(line[30:38])
y = float(line[38:46])
z = float(line[46:54])
x_max = max(x, x_max)
y_max = max(y, y_max)
z_max = max(z, z_max)
x_min = min(x, x_min)
y_min = min(y, y_min)
z_min = min(z, z_min)
center_x = (x_min + x_max) / 2.0
center_y = (y_min + y_max) / 2.0
center_z = (z_min + z_max) / 2.0
size_x = x_max - x_min + 2 * padding
size_y = y_max - y_min + 2 * padding
size_z = z_max - z_min + 2 * padding
return (center_x, center_y, center_z), (size_x, size_y, size_z)

boron_silicon_atompar = "atom_par Si 4.10 0.200 35.8235 -0.00143 0.0 0.0 0 -1 -1 6" + os_linesep
boron_silicon_atompar += "atom_par B 3.84 0.155 29.6478 -0.00152 0.0 0.0 0 -1 -1 0" + os_linesep
33 changes: 17 additions & 16 deletions meeko/receptor_pdbqt.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
import json
from os import linesep as os_linesep
import pathlib
import sys

import numpy as np
from scipy import spatial
Expand Down Expand Up @@ -86,7 +87,7 @@ def _read_receptor_pdbqt_string(pdbqt_string, skip_typing=False):
except:
temp_factor = None
record_type = line[0:6].strip()

if skip_typing:
atoms.append((idx, serial, name, resid, resname, chainid, xyz, partial_charges, atom_type,
alt_id, in_code, occupancy, temp_factor, record_type))
Expand Down Expand Up @@ -135,7 +136,7 @@ def __init__(self, pdbqt_filename, skip_typing=False):
self._KDTree = None

with open(pdbqt_filename) as f:
pdbqt_string = f.read()
pdbqt_string = f.read()

self._atoms, self._atom_annotations = _read_receptor_pdbqt_string(pdbqt_string, skip_typing)
# We add to the KDTree only the rigid part of the receptor
Expand All @@ -150,7 +151,7 @@ def __repr__(self):
def get_atom_indices_by_residue(atoms):
""" return a dictionary where residues are keys and
values are lists of atom indices
>>> atom_idx_by_res = {("A", "LYS", 417): [0, 1, 2, 3, ..., 8]}
"""

Expand Down Expand Up @@ -184,18 +185,18 @@ def get_params_for_residue(resname, atom_names, residue_params=residue_params):
if not is_matched:
ok = False
return atom_params, ok, err

for atom_name in atom_names:
name_index = residue_params[r_id]["atom_names"].index(atom_name)
for param in residue_params[r_id].keys():
if param in excluded_params:
continue
if param not in atom_params:
atom_params[param] = [None] * atom_counter
atom_params[param] = [None] * atom_counter
value = residue_params[r_id][param][name_index]
atom_params[param].append(value)
atom_counter += 1

return atom_params, ok, err

def assign_types_charges(self, residue_params=residue_params):
Expand All @@ -215,8 +216,8 @@ def assign_types_charges(self, residue_params=residue_params):
for key in wanted_params:
atom_params[key].extend(params_this_res[key])
if ok:
self._atoms["partial_charges"] = atom_params["gasteiger"]
self._atoms["atom_type"] = atom_params["atom_types"]
self._atoms["partial_charges"] = atom_params["gasteiger"]
self._atoms["atom_type"] = atom_params["atom_types"]
return ok, err

def write_flexres_from_template(self, res_id, atom_index=0):
Expand Down Expand Up @@ -259,7 +260,7 @@ def write_flexres_from_template(self, res_id, atom_index=0):
for i in range(len(template["is_atom"])):
if template["is_atom"][i]:
ref_atoms.add(template["atom_name"][i])
if got_atoms != ref_atoms:
if got_atoms != ref_atoms:
success = False
error_msg += "mismatch in atom names for residue %s" % str(res_id) + os_linesep
error_msg += "names found but not in template: %s" % str(got_atoms.difference(ref_atoms)) + os_linesep
Expand All @@ -273,7 +274,7 @@ def write_flexres_from_template(self, res_id, atom_index=0):
atom_index += 1
name = template['atom_name'][i]
atom = atoms_by_name[name]
if atom["atom_type"] not in self.skip_types:
if atom["atom_type"] not in self.skip_types:
atom["serial"] = atom_index
output["pdbqt"] += self.write_pdbqt_line(atom)
else:
Expand Down Expand Up @@ -319,7 +320,7 @@ def write_pdbqt_string(self, flexres=()):
for i, atom in enumerate(self._atoms):
if i not in pdbqt["flex_indices"] and atom["atom_type"] not in self.skip_types:
pdbqt["rigid"] += self.write_pdbqt_line(atom)

return pdbqt, ok, err

@staticmethod
Expand All @@ -343,7 +344,7 @@ def get_neigh(idx, bonds):
two_bond_away.add(j)
names_1bond = [atom_names[i] for i in one_bond_away]
names_2bond = [atom_names[i] for i in two_bond_away]
new_pdbqt_str = ""
new_pdbqt_str = ""
for i, line in enumerate(pdbqtstr.split(os_linesep)[:-1]):
if line.startswith("ATOM") or line.startswith("HETATM"):
name = line[12:16].strip()
Expand Down Expand Up @@ -394,13 +395,13 @@ def positions(self, atom_idx=None):
return np.atleast_2d(self.atoms(atom_idx)['xyz'])

def closest_atoms_from_positions(self, xyz, radius, atom_properties=None, ignore=None):
"""Retrieve indices of the closest atoms around a positions/coordinates
"""Retrieve indices of the closest atoms around a positions/coordinates
at a certain radius.
Args:
xyz (np.ndarray): array of 3D coordinates
raidus (float): radius
atom_properties (str): property of the atoms to retrieve
atom_properties (str): property of the atoms to retrieve
(properties: ligand, flexible_residue, vdw, hb_don, hb_acc, metal, water, reactive, glue)
ignore (int or list): ignore atom for the search using atom id (0-based)
Expand Down Expand Up @@ -442,13 +443,13 @@ def closest_atoms_from_positions(self, xyz, radius, atom_properties=None, ignore
return atoms

def closest_atoms(self, atom_idx, radius, atom_properties=None):
"""Retrieve indices of the closest atoms around a positions/coordinates
"""Retrieve indices of the closest atoms around a positions/coordinates
at a certain radius.
Args:
atom_idx (int, list): index of one or multiple atoms (0-based)
raidus (float): radius
atom_properties (str or list): property of the atoms to retrieve
atom_properties (str or list): property of the atoms to retrieve
(properties: ligand, flexible_residue, vdw, hb_don, hb_acc, metal, water, reactive, glue)
Returns:
Expand Down
67 changes: 40 additions & 27 deletions scripts/mk_prepare_receptor.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ def parse_residue_string(string):
err += "resnum could not be converted to integer, it was '%s' in '%s'" % (resnum, string) + os_linesep
if ok:
res_id = (chain, resname, resnum)
return res_id, ok, err
return res_id, ok, err

def parse_residue_string_and_name(string):
"""
Expand Down Expand Up @@ -127,6 +127,8 @@ def get_args():
parser.add_argument('--box_size', help="size of grid box (x, y, z) in Angstrom", nargs=3, type=float)
parser.add_argument('--box_center', help="center of grid box (x, y, z) in Angstrom", nargs=3, type=float)
parser.add_argument('--box_center_on_reactive_res', help="project center of grid box along CA-CB bond 5 A away from CB", action="store_true")
parser.add_argument('--ligand', help="PDBQT of reference ligand")
parser.add_argument('--padding', help="padding around reference ligand [A]", type=float)
parser.add_argument('--skip_gpf', help="do not write a GPF file for autogrid", action="store_true")
parser.add_argument('--r_eq_12', default=1.8, type=float, help="r_eq for reactive atoms (1-2 interaction)")
parser.add_argument('--eps_12', default=2.5, type=float, help="epsilon for reactive atoms (1-2 interaction)")
Expand All @@ -142,21 +144,28 @@ def get_args():
msg = "can't use both --box_center and --box_center_on_reactive_res"
print("Command line error: " + msg, file=sys.stderr)
sys.exit(2)
got_center = (args.box_center is not None) or args.box_center_on_reactive_res
if not args.skip_gpf and (args.box_size is None) == got_center:
msg = "missing center or size of grid box to write .gpf file for autogrid4" + os_linesep
msg += "use --box_size and either --box_center or --box_center_on_reactive_res" + os_linesep
msg += "Exactly one reactive residue required for --box_center_on_reactive_res" + os_linesep
msg += "If a GPF file is not needed (e.g. docking with Vina scoring function) use option --skip_gpf"
print("Command line error: " + msg, file=sys.stderr)
sys.exit(2)
if (args.box_size is None) and not args.skip_gpf:
msg = "grid box information is needed to dock with the AD4 scoring function." + os_linesep
msg += "The grid box center and size will be used to write a GPF file for autogrid" + os_linesep
msg += "If a GPF file is not needed (e.g. docking with Vina scoring function) use option --skip_gpf"
got_center = (args.box_center is not None) or args.box_center_on_reactive_res or (args.ligand is not None)
if not args.skip_gpf:
if not got_center:
msg = "missing center or size of grid box to write .gpf file for autogrid4" + os_linesep
msg += "use --box_size and either --box_center or --box_center_on_reactive_res" + os_linesep
msg += "or --ligand and --padding" + os_linesep
msg += "Exactly one reactive residue required for --box_center_on_reactive_res" + os_linesep
msg += "If a GPF file is not needed (e.g. docking with Vina scoring function) use option --skip_gpf"
print("Command line error: " + msg, file=sys.stderr)
sys.exit(2)
if (args.box_size is None) and (args.padding is None):
msg = "grid box information is needed to dock with the AD4 scoring function." + os_linesep
msg += "The grid box center and size will be used to write a GPF file for autogrid" + os_linesep
msg += "If a GPF file is not needed (e.g. docking with Vina scoring function) use option --skip_gpf"
print("Command line error: " + msg, file=sys.stderr)
sys.exit(2)

if (args.box_center is not None) + (args.ligand is not None) + args.box_center_on_reactive_res > 1:
msg = "--box_center, --box_center_on_reactive_res, and --ligand are mutually exclusive options"
print("Command line error: " + msg, file=sys.stderr)
sys.exit(2)

return args

args = get_args()
Expand All @@ -183,34 +192,34 @@ def get_args():
all_ok = True
all_err = ""
for resid_string in args.reactive_flexres:
res_id, ok, err = parse_residue_string(resid_string)
res_id, ok, err = parse_residue_string(resid_string)
if ok:
resname = res_id[1]
if resname in reactive_atom:
reactive_flexres[res_id] = reactive_atom[resname]
else:
all_ok = False
all_err += "no default reactive name for %s, " % resname
all_err += "no default reactive name for %s, " % resname
all_err += "use --reactive_name or --reactive_name_specific" + os_linesep
all_ok &= ok
all_err += err

for string in args.reactive_name_specific:
out, ok, err = parse_residue_string_and_name(string)
out, ok, err = parse_residue_string_and_name(string)
if ok:
# override name if res_id was also passed to --reactive_flexres
reactive_flexres[out["res_id"]] = out["name"]
all_ok &= ok
all_err += err

if len(reactive_flexres) > 8:
if len(reactive_flexres) > 8:
msg = "got %d reactive_flexres but maximum is 8." % (len(args.reactive_flexres))
print("Command line error: " + msg, file=sys.stderr)
sys.exit(2)

all_flexres = set()
for resid_string in args.flexres:
res_id, ok, err = parse_residue_string(resid_string)
res_id, ok, err = parse_residue_string(resid_string)
all_ok &= ok
all_err += err
if ok:
Expand Down Expand Up @@ -247,7 +256,7 @@ def get_args():
msg += "residue, but %d reactive residues are set" % len(reactive_flexres)
print("Command line error:" + msg, file=sys.stderr)
sys.exit(2)

if args.pdb is not None:
receptor = PDBQTReceptor(args.pdb, skip_typing=True)
ok, err = receptor.assign_types_charges()
Expand All @@ -257,7 +266,7 @@ def get_args():

any_lig_base_types = ["HD", "C", "A", "N", "NA", "OA", "F", "P", "SA",
"S", "Cl", "CL", "Br", "BR", "I", "Si", "B"]

pdbqt, ok, err = receptor.write_pdbqt_string(flexres=all_flexres)
check(ok, err)

Expand All @@ -283,7 +292,7 @@ def get_args():

suffix = outpath.suffix
if outpath.suffix == "":
suffix = ".pdbqt"
suffix = ".pdbqt"
rigid_fn = str(outpath.with_suffix("")) + "_rigid" + suffix
flex_fn = str(outpath.with_suffix("")) + "_flex" + suffix

Expand All @@ -295,16 +304,17 @@ def get_args():
written_files_log["filename"].append(rigid_fn)
written_files_log["description"].append("static (i.e., rigid) receptor input file")
with open(rigid_fn, "w") as f:
f.write(pdbqt["rigid"])
f.write(pdbqt["rigid"])

# GPF for autogrid4
if not args.skip_gpf:
if args.box_center is not None:
box_center = args.box_center
box_size = args.box_size
elif args.box_center_on_reactive_res:
# we have only one reactive residue and will set the box center
# to be 5 Angstromg away from CB along the CA->CB vector
idxs = receptor.atom_idxs_by_res[list(reactive_flexres.keys())[0]]
idxs = receptor.atom_idxs_by_res[list(reactive_flexres.keys())[0]]
ca = None
cb = None
for atom in receptor.atoms(idxs):
Expand All @@ -317,6 +327,9 @@ def get_args():
v = (cb - ca)
v /= math.sqrt(v[0]**2 + v[1]**2 + v[2]**2) + 1e-8
box_center = ca + 5 * v
box_size = args.box_size
elif args.ligand is not None:
box_center, box_size = gridbox.calc_box(args.ligand, args.padding)
else:
print("Error: No box center specified.", file=sys.stderr)
sys.exit(2)
Expand All @@ -328,7 +341,7 @@ def get_args():
with open(ff_fn, "w") as f:
f.write(gridbox.boron_silicon_atompar)
rec_types = set(t for (i, t) in enumerate(receptor.atoms()["atom_type"]) if i not in pdbqt["flex_indices"])
gpf_string, npts = gridbox.get_gpf_string(box_center, args.box_size, rigid_fn, rec_types, any_lig_base_types,
gpf_string, npts = gridbox.get_gpf_string(box_center, box_size, rigid_fn, rec_types, any_lig_base_types,
ff_param_fname=ff_fn.name)
# write GPF
gpf_fn = pathlib.Path(rigid_fn).with_suffix(".gpf")
Expand All @@ -351,7 +364,7 @@ def get_args():
print("WARNING: Flexible residue outside box." + os_linesep, file=sys.stderr)
print("WARNING: Strongly recommended to use a box that encompasses flexible residues." + os_linesep, file=sys.stderr)
break # only need to warn once

# configuration info for AutoDock-GPU reactive docking
if len(reactive_flexres) > 0:
any_lig_reac_types = []
Expand All @@ -363,7 +376,7 @@ def get_args():
for line in all_flex_pdbqt.split(os_linesep):
if line.startswith("ATOM") or line.startswith("HETATM"):
atype = line[77:].strip()
basetype, _ = reactive_typer.get_basetype_and_order(atype)
basetype, _ = reactive_typer.get_basetype_and_order(atype)
if basetype is not None: # is None if not reactive
rec_reac_types.append(line[77:].strip())

Expand Down

0 comments on commit efb14bc

Please sign in to comment.