OpenStructure
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
dockq.py
Go to the documentation of this file.
1 import sys
2 import os
3 import subprocess
4 import tempfile
5 import shutil
6 
7 from ost import io
8 from ost import mol
9 
10 def _Setup(mdl, ref, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2):
11  """ Performs parameter checks and dumps files for DockQ
12 
13  In case of dimeric interfaces the respective chains are selected from
14  mdl/trg, renamed to A/B and dumped to disk.
15 
16  In case of interfaces with more chains involved, we simply select the
17  specified chains and do no renaming before dumping to disk.
18  """
19  if isinstance(mdl_ch1, str):
20  mdl_ch1 = [mdl_ch1]
21  if isinstance(mdl_ch2, str):
22  mdl_ch2 = [mdl_ch2]
23  if isinstance(ref_ch1, str):
24  ref_ch1 = [ref_ch1]
25  if isinstance(ref_ch2, str):
26  ref_ch2 = [ref_ch2]
27 
28  if len(mdl_ch1) == 0:
29  raise RuntimeError("mdl_ch1 is empty")
30  if len(mdl_ch2) == 0:
31  raise RuntimeError("mdl_ch2 is empty")
32 
33  if len(mdl_ch1) != len(ref_ch1):
34  raise RuntimeError("mdl_ch1/ref_ch1 inconsistent in size")
35  if len(mdl_ch2) != len(ref_ch2):
36  raise RuntimeError("mdl_ch2/ref_ch2 inconsistent in size")
37 
38  for cname in mdl_ch1:
39  ch = mdl.FindChain(cname)
40  if not ch.IsValid():
41  raise RuntimeError(f"Chain {cname} specified in mdl_ch1 not "
42  f"present in mdl")
43 
44  for cname in mdl_ch2:
45  ch = mdl.FindChain(cname)
46  if not ch.IsValid():
47  raise RuntimeError(f"Chain {cname} specified in mdl_ch2 not "
48  f"present in mdl")
49 
50  for cname in ref_ch1:
51  ch = ref.FindChain(cname)
52  if not ch.IsValid():
53  raise RuntimeError(f"Chain {cname} specified in ref_ch1 not "
54  f"present in ref")
55 
56  for cname in ref_ch2:
57  ch = ref.FindChain(cname)
58  if not ch.IsValid():
59  raise RuntimeError(f"Chain {cname} specified in ref_ch2 not "
60  f"present in ref")
61 
62  mdl_to_dump = mdl.CreateFullView()
63  ref_to_dump = ref.CreateFullView()
64 
65  if len(mdl_ch1) == 1 and len(mdl_ch2) == 1:
66  # Dimer processing of mdl => Create new entity only containing
67  # the two specified chains and rename them to A, B
68  mdl_to_dump = mol.CreateEntityFromView(mdl_to_dump, True)
69  tmp = mol.CreateEntity()
70  ed = tmp.EditXCS()
71  ch1 = mdl_to_dump.FindChain(mdl_ch1[0])
72  ed.InsertChain("A", ch1, deep=True)
73  ch2 = mdl_to_dump.FindChain(mdl_ch2[0])
74  ed.InsertChain("B", ch2, deep=True)
75  mdl_ch1 = ["A"]
76  mdl_ch2 = ["B"]
77  mdl_to_dump = tmp
78 
79  # Same for ref
80  ref_to_dump = mol.CreateEntityFromView(ref_to_dump, True)
81  tmp = mol.CreateEntity()
82  ed = tmp.EditXCS()
83  ch1 = ref_to_dump.FindChain(ref_ch1[0])
84  ed.InsertChain("A", ch1, deep=True)
85  ch2 = ref_to_dump.FindChain(ref_ch2[0])
86  ed.InsertChain("B", ch2, deep=True)
87  ref_ch1 = ["A"]
88  ref_ch2 = ["B"]
89  ref_to_dump = tmp
90  else:
91  # Interface with more chains...
92  raise NotImplementedError("DockQ computations beyond two interacting "
93  "chains has not been properly tested...")
94  #mdl_chain_names = mdl_ch1 + mdl_ch2
95  #ref_chain_names = ref_ch1 + ref_ch2
96  #mdl_to_dump = mdl_to_dump.Select(f"cname={','.join(mdl_chain_names)}")
97  #ref_to_dump = ref_to_dump.Select(f"cname={','.join(ref_chain_names)}")
98 
99  # first write structures to string, only create a tmpdir and the actual
100  # files if this succeeds
101  mdl_str = io.EntityToPDBStr(mdl_to_dump)
102  ref_str = io.EntityToPDBStr(ref_to_dump)
103 
104  tmp_dir = tempfile.mkdtemp()
105  with open(os.path.join(tmp_dir, "mdl.pdb"), 'w') as fh:
106  fh.write(mdl_str)
107  with open(os.path.join(tmp_dir, "ref.pdb"), 'w') as fh:
108  fh.write(ref_str)
109 
110  return (tmp_dir, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2)
111 
113  """ DockQ result object
114  """
115  def __init__(self, Fnat, Fnonnat, native_contacts, model_contacts, iRMS,
116  LRMS, DockQ):
117  self._Fnat = Fnat
118  self._Fnonnat = Fnonnat
119  self._native_contacts = native_contacts
120  self._model_contacts = model_contacts
121  self._iRMS = iRMS
122  self._LRMS = LRMS
123  self._DockQ = DockQ
124 
125  @property
126  def Fnat(self):
127  """ DockQ - Fnat output
128 
129  :type: :class:`float`
130  """
131  return self._Fnat
132 
133  @property
134  def Fnonnat(self):
135  """ DockQ - Fnonnat output
136 
137  :type: :class:`float`
138  """
139  return self._Fnonnat
140 
141  @property
142  def native_contacts(self):
143  """ DockQ - number native contacts
144 
145  :type: :class:`int`
146  """
147  return self._native_contacts
148 
149  @property
150  def model_contacts(self):
151  """ DockQ - number model contacts
152 
153  :type: :class:`int`
154  """
155  return self._model_contacts
156 
157  @property
158  def iRMS(self):
159  """ DockQ - iRMS output
160 
161  :type: :class:`float`
162  """
163  return self._iRMS
164 
165  @property
166  def LRMS(self):
167  """ DockQ - LMRS output
168 
169  :type: :class:`float`
170  """
171  return self._LRMS
172 
173  @property
174  def DockQ(self):
175  """ DockQ - DockQ output
176 
177  :type: :class:`float`
178  """
179  return self._DockQ
180 
181  def JSONSummary(self):
182  """ Returns JSON serializable summary
183  """
184  return {"Fnat": self.Fnat,
185  "Fnonnat": self.Fnonnat,
186  "native_contacts": self.native_contacts,
187  "model_contacts": self.model_contacts,
188  "iRMS": self.iRMS,
189  "LRMS": self.LRMS,
190  "DockQ": self.DockQ}
191 
192  @staticmethod
193  def FromDockQOutput(output):
194  """ Static constructor from raw DockQ output
195 
196  :param output: Raw output from DockQ executable
197  :type output: :class:`str`
198  :returns: Object of type :class:`DockQResult`
199  """
200  Fnat = None
201  Fnonnat = None
202  native_contacts = None
203  model_contacts = None
204  iRMS = None
205  LRMS = None
206  DockQ = None
207 
208  for line in output.splitlines():
209  if line.startswith('*'):
210  continue
211  if line.startswith("Fnat"):
212  Fnat = float(line.split()[1])
213  native_contacts = int(line.split()[5])
214  elif line.startswith("Fnonnat"):
215  Fnonnat = float(line.split()[1])
216  model_contacts = int(line.split()[5])
217  elif line.startswith("iRMS"):
218  iRMS = float(line.split()[1])
219  elif line.startswith("LRMS"):
220  LRMS = float(line.split()[1])
221  elif line.startswith("DockQ"):
222  DockQ = float(line.split()[1])
223 
224  return DockQResult(Fnat, Fnonnat, native_contacts, model_contacts,
225  iRMS, LRMS, DockQ)
226 
227 
228 def DockQ(dockq_exec, mdl, ref, mdl_ch1, mdl_ch2, ref_ch1,
229  ref_ch2):
230  """ Computes DockQ for specified interface
231 
232  DockQ is available from https://github.com/bjornwallner/DockQ -
233  For this binding to work, DockQ must be properly installed and its
234  dependencies must be available (numpy, Biopython).
235 
236  :param dockq_exec: Path to DockQ.py script from DockQ repository
237  :type dockq_exec: :class:`str`
238  :param mdl: Model structure
239  :type mdl: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
240  :param ref: Reference structure, i.e. native structure
241  :type ref: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
242  :param mdl_ch1: Specifies chain(s) in model constituting first part of
243  interface
244  :type mdl_ch1: :class:`str`/:class:`list` of :class:`str`
245  :param mdl_ch2: Specifies chain(s) in model constituting second part of
246  interface
247  :type mdl_ch2: :class:`str`/:class:`list` of :class:`str`
248  :param ref_ch1: ref equivalent of mdl_ch1
249  :type ref_ch1: :class:`str`/:class:`list` of :class:`str`
250  :param ref_ch2: ref equivalent of mdl_ch2
251  :type ref_ch2: :class:`str`/:class:`list` of :class:`str`
252  :returns: Result object of type :class:`DockQResult`
253  """
254  if not os.path.exists(dockq_exec):
255  raise RuntimeError(f"DockQ executable ({dockq_exec}) does not exist")
256 
257  tmp_dir, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2 = \
258  _Setup(mdl, ref, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2)
259 
260  cmd = [sys.executable, dockq_exec, os.path.join(tmp_dir, "mdl.pdb"),
261  os.path.join(tmp_dir, "ref.pdb")]
262 
263  # add mdl/ref chains
264  cmd.append("-model_chain1")
265  cmd += mdl_ch1
266  cmd.append("-model_chain2")
267  cmd += mdl_ch2
268  cmd.append("-native_chain1")
269  cmd += ref_ch1
270  cmd.append("-native_chain2")
271  cmd += ref_ch2
272 
273  proc = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
274 
275  shutil.rmtree(tmp_dir) # cleanup, no matter if DockQ executed successfully
276 
277  if proc.returncode != 0:
278  raise RuntimeError("DockQ run failed - returncode: " + \
279  str(proc.returncode) + ", stderr: " + \
280  proc.stderr.decode() + ", stdout: " + \
281  proc.stdout.decode())
282 
283  if proc.stderr.decode() != "":
284  raise RuntimeError("DockQ run failed - stderr: " + \
285  proc.stderr.decode() + ", stdout: " + \
286  proc.stdout.decode())
287 
288  return DockQResult.FromDockQOutput(proc.stdout.decode())