File: dockq.py

package info (click to toggle)
openstructure 2.11.1-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 206,240 kB
  • sloc: cpp: 188,571; python: 36,686; ansic: 34,298; fortran: 3,275; sh: 312; xml: 146; makefile: 29
file content (284 lines) | stat: -rw-r--r-- 9,258 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
import sys
import os
import subprocess
import tempfile
import shutil

from ost import io
from ost import mol

def _Setup(mdl, ref, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2):
    """ Performs parameter checks and dumps files for DockQ

    In case of dimeric interfaces the respective chains are selected from
    mdl/trg, renamed to A/B and dumped to disk.

    In case of interfaces with more chains involved, we simply select the
    specified chains and do no renaming before dumping to disk.
    """
    if isinstance(mdl_ch1, str):
        mdl_ch1 = [mdl_ch1]
    if isinstance(mdl_ch2, str):
        mdl_ch2 = [mdl_ch2]
    if isinstance(ref_ch1, str):
        ref_ch1 = [ref_ch1]
    if isinstance(ref_ch2, str):
        ref_ch2 = [ref_ch2]

    if len(mdl_ch1) == 0:
        raise RuntimeError("mdl_ch1 is empty")
    if len(mdl_ch2) == 0:
        raise RuntimeError("mdl_ch2 is empty")

    if len(mdl_ch1) != len(ref_ch1):
        raise RuntimeError("mdl_ch1/ref_ch1 inconsistent in size")
    if len(mdl_ch2) != len(ref_ch2):
        raise RuntimeError("mdl_ch2/ref_ch2 inconsistent in size")

    for cname in mdl_ch1:
        ch = mdl.FindChain(cname)
        if not ch.IsValid():
            raise RuntimeError(f"Chain {cname} specified in mdl_ch1 not "
                               f"present in mdl")

    for cname in mdl_ch2:
        ch = mdl.FindChain(cname)
        if not ch.IsValid():
            raise RuntimeError(f"Chain {cname} specified in mdl_ch2 not "
                               f"present in mdl")

    for cname in ref_ch1:
        ch = ref.FindChain(cname)
        if not ch.IsValid():
            raise RuntimeError(f"Chain {cname} specified in ref_ch1 not "
                               f"present in ref")

    for cname in ref_ch2:
        ch = ref.FindChain(cname)
        if not ch.IsValid():
            raise RuntimeError(f"Chain {cname} specified in ref_ch2 not "
                               f"present in ref")

    mdl_to_dump = mdl.CreateFullView()
    ref_to_dump = ref.CreateFullView()

    if len(mdl_ch1) == 1 and len(mdl_ch2) == 1:
        # Dimer processing of mdl => Create new entity only containing 
        # the two specified chains and rename them to A, B
        mdl_to_dump = mol.CreateEntityFromView(mdl_to_dump, True)
        tmp = mol.CreateEntity()
        ed = tmp.EditXCS()
        ch1 = mdl_to_dump.FindChain(mdl_ch1[0])
        ed.InsertChain("A", ch1, deep=True)
        ch2 = mdl_to_dump.FindChain(mdl_ch2[0])
        ed.InsertChain("B", ch2, deep=True)
        mdl_ch1 = ["A"]
        mdl_ch2 = ["B"]
        mdl_to_dump = tmp

        # Same for ref
        ref_to_dump = mol.CreateEntityFromView(ref_to_dump, True)
        tmp = mol.CreateEntity()
        ed = tmp.EditXCS()
        ch1 = ref_to_dump.FindChain(ref_ch1[0])
        ed.InsertChain("A", ch1, deep=True)
        ch2 = ref_to_dump.FindChain(ref_ch2[0])
        ed.InsertChain("B", ch2, deep=True)
        ref_ch1 = ["A"]
        ref_ch2 = ["B"]
        ref_to_dump = tmp
    else:
        # Interface with more chains...
        raise NotImplementedError("DockQ computations beyond two interacting "
                                  "chains has not been properly tested...")

    # first write structures to string, only create a tmpdir and the actual
    # files if this succeeds
    mdl_str = io.EntityToPDBStr(mdl_to_dump)
    ref_str = io.EntityToPDBStr(ref_to_dump)

    tmp_dir = tempfile.mkdtemp()
    with open(os.path.join(tmp_dir, "mdl.pdb"), 'w') as fh:
        fh.write(mdl_str)
    with open(os.path.join(tmp_dir, "ref.pdb"), 'w') as fh:
        fh.write(ref_str)

    return (tmp_dir, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2)

class DockQResult:
    """ DockQ result object
    """
    def __init__(self, Fnat, Fnonnat, native_contacts, model_contacts, iRMS,
                 LRMS, DockQ):
        self._Fnat = Fnat
        self._Fnonnat = Fnonnat
        self._native_contacts = native_contacts
        self._model_contacts = model_contacts
        self._iRMS = iRMS
        self._LRMS = LRMS
        self._DockQ = DockQ

    @property
    def Fnat(self):
        """ DockQ - Fnat output

        :type: :class:`float`
        """
        return self._Fnat
    
    @property
    def Fnonnat(self):
        """ DockQ - Fnonnat output

        :type: :class:`float`
        """
        return self._Fnonnat

    @property
    def native_contacts(self):
        """ DockQ - number native contacts

        :type: :class:`int`
        """
        return self._native_contacts

    @property
    def model_contacts(self):
        """ DockQ - number model contacts

        :type: :class:`int`
        """
        return self._model_contacts

    @property
    def iRMS(self):
        """ DockQ - iRMS output

        :type: :class:`float`
        """
        return self._iRMS

    @property
    def LRMS(self):
        """ DockQ - LMRS output 

        :type: :class:`float`
        """
        return self._LRMS

    @property
    def DockQ(self):
        """ DockQ - DockQ output

        :type: :class:`float`
        """
        return self._DockQ

    def JSONSummary(self):
        """ Returns JSON serializable summary
        """
        return {"Fnat": self.Fnat,
                "Fnonnat": self.Fnonnat,
                "native_contacts": self.native_contacts,
                "model_contacts": self.model_contacts,
                "iRMS": self.iRMS,
                "LRMS": self.LRMS,
                "DockQ": self.DockQ}

    @staticmethod
    def FromDockQOutput(output):
        """ Static constructor from raw DockQ output

        :param output: Raw output from DockQ executable
        :type output: :class:`str`
        :returns: Object of type :class:`DockQResult`
        """
        Fnat = None
        Fnonnat = None
        native_contacts = None
        model_contacts = None
        iRMS = None
        LRMS = None
        DockQ = None

        for line in output.splitlines():
            if line.startswith('*'):
                continue
            if line.startswith("Fnat"):
                Fnat = float(line.split()[1])
                native_contacts = int(line.split()[5])
            elif line.startswith("Fnonnat"):
                Fnonnat = float(line.split()[1])
                model_contacts = int(line.split()[5])
            elif line.startswith("iRMS"):
                iRMS = float(line.split()[1])
            elif line.startswith("LRMS"):
                LRMS = float(line.split()[1])
            elif line.startswith("DockQ"):
                DockQ = float(line.split()[1])

        return DockQResult(Fnat, Fnonnat, native_contacts, model_contacts,
                           iRMS, LRMS, DockQ)


def DockQ(dockq_exec, mdl, ref, mdl_ch1, mdl_ch2, ref_ch1,
          ref_ch2):
    """ Computes DockQ for specified interface 

    DockQ is available from https://github.com/bjornwallner/DockQ - 
    For this binding to work, DockQ must be properly installed and its
    dependencies must be available (numpy, Biopython).

    :param dockq_exec: Path to DockQ.py script from DockQ repository
    :type dockq_exec: :class:`str`
    :param mdl: Model structure
    :type mdl: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
    :param ref: Reference structure, i.e. native structure
    :type ref: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
    :param mdl_ch1: Specifies chain(s) in model constituting first part of
                    interface
    :type mdl_ch1: :class:`str`/:class:`list` of :class:`str`
    :param mdl_ch2: Specifies chain(s) in model constituting second part of
                    interface
    :type mdl_ch2: :class:`str`/:class:`list` of :class:`str`
    :param ref_ch1: ref equivalent of mdl_ch1
    :type ref_ch1: :class:`str`/:class:`list` of :class:`str`
    :param ref_ch2: ref equivalent of mdl_ch2
    :type ref_ch2: :class:`str`/:class:`list` of :class:`str`
    :returns: Result object of type :class:`DockQResult`
    """
    if not os.path.exists(dockq_exec):
        raise RuntimeError(f"DockQ executable ({dockq_exec}) does not exist")

    tmp_dir, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2 = \
    _Setup(mdl, ref, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2)

    cmd = [sys.executable, dockq_exec, os.path.join(tmp_dir, "mdl.pdb"),
           os.path.join(tmp_dir, "ref.pdb")]

    # add mdl/ref chains
    cmd.append("-model_chain1")
    cmd += mdl_ch1
    cmd.append("-model_chain2")
    cmd += mdl_ch2
    cmd.append("-native_chain1")
    cmd += ref_ch1
    cmd.append("-native_chain2")
    cmd += ref_ch2

    proc = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)

    shutil.rmtree(tmp_dir) # cleanup, no matter if DockQ executed successfully

    if proc.returncode != 0:
        raise RuntimeError("DockQ run failed - returncode: " + \
                           str(proc.returncode) + ", stderr: " + \
                           proc.stderr.decode() + ", stdout: " + \
                           proc.stdout.decode())

    if proc.stderr.decode() != "":
        raise RuntimeError("DockQ run failed - stderr: " + \
                           proc.stderr.decode() + ", stdout: " + \
                           proc.stdout.decode())

    return DockQResult.FromDockQOutput(proc.stdout.decode())