diff --git a/modules/mol/alg/pymod/dockq.py b/modules/mol/alg/pymod/dockq.py index c22ca78d21746de90900b2a1f539d931c8cc65a8..126c8f57519472b8cb74ccaa45e5b0719c9043f8 100644 --- a/modules/mol/alg/pymod/dockq.py +++ b/modules/mol/alg/pymod/dockq.py @@ -1,5 +1,6 @@ from ost import geom from ost import mol +from ost import seq def _PreprocessStructures(mdl, ref, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2, ch1_aln = None, ch2_aln = None): @@ -15,7 +16,56 @@ def _PreprocessStructures(mdl, ref, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2, ref_residues_1 = list() ref_residues_2 = list() - if ch1_aln is None and ch2_aln is None: + if ch1_aln is not None and ch2_aln is not None: + # there are potentially already views attached to the alns but + # we go for *mdl* and *ref* here + if ch1_aln.GetCount() != 2 or ch2_aln.GetCount() != 2: + raise RuntimeError("Expect exactly two sequences in provided alns!") + + tmp = ch1_aln.GetSequence(0) + ref_s1 = seq.CreateSequence(tmp.GetName(), str(tmp)) + ref_s1.SetOffset(tmp.GetOffset()) + ref_s1.AttachView(ref.Select(f"cname={ref_ch1}")) + tmp = ch1_aln.GetSequence(1) + mdl_s1 = seq.CreateSequence(tmp.GetName(), str(tmp)) + mdl_s1.SetOffset(tmp.GetOffset()) + mdl_s1.AttachView(mdl.Select(f"cname={mdl_ch1}")) + new_ch1_aln = seq.CreateAlignment(ref_s1, mdl_s1) + for col in new_ch1_aln: + if col[0] != '-' and col[1] != '-': + ref_r = col.GetResidue(0) + mdl_r = col.GetResidue(1) + if not (ref_r.IsValid() and ref_r.one_letter_code == col[0]): + raise RuntimeError("DockQ: mismatch between provided " + "alignments and ATOMSEQ in structures") + if not (mdl_r.IsValid() and mdl_r.one_letter_code == col[1]): + raise RuntimeError("DockQ: mismatch between provided " + "alignments and ATOMSEQ in structures") + ref_residues_1.append(ref_r) + mdl_residues_1.append(mdl_r) + + tmp = ch2_aln.GetSequence(0) + ref_s2 = seq.CreateSequence(tmp.GetName(), str(tmp)) + ref_s2.SetOffset(tmp.GetOffset()) + ref_s2.AttachView(ref.Select(f"cname={ref_ch2}")) + tmp = ch2_aln.GetSequence(1) + mdl_s2 = seq.CreateSequence(tmp.GetName(), str(tmp)) + mdl_s2.SetOffset(tmp.GetOffset()) + mdl_s2.AttachView(mdl.Select(f"cname={mdl_ch2}")) + new_ch2_aln = seq.CreateAlignment(ref_s2, mdl_s2) + for col in new_ch2_aln: + if col[0] != '-' and col[1] != '-': + ref_r = col.GetResidue(0) + mdl_r = col.GetResidue(1) + if not (ref_r.IsValid() and ref_r.one_letter_code == col[0]): + raise RuntimeError("DockQ: mismatch between provided " + "alignments and ATOMSEQ in structures") + if not (mdl_r.IsValid() and mdl_r.one_letter_code == col[1]): + raise RuntimeError("DockQ: mismatch between provided " + "alignments and ATOMSEQ in structures") + ref_residues_2.append(ref_r) + mdl_residues_2.append(mdl_r) + else: # go by residue numbers for mdl_r in mdl.Select(f"cname={mdl_ch1}").residues: ref_r = ref.FindResidue(ref_ch1, mdl_r.GetNumber()) @@ -27,8 +77,6 @@ def _PreprocessStructures(mdl, ref, mdl_ch1, mdl_ch2, ref_ch1, ref_ch2, if ref_r.IsValid(): mdl_residues_2.append(mdl_r) ref_residues_2.append(ref_r) - else: - raise NotImplementedError("No aln mapping implemented yet") new_mdl = mdl.handle.CreateEmptyView() new_ref = ref.handle.CreateEmptyView() diff --git a/modules/mol/alg/pymod/scoring.py b/modules/mol/alg/pymod/scoring.py index 8c484570e612dd07a706efe15452ae0a4efdbef2..5555282b593083a3e3bc90282ead1db01bc52d1f 100644 --- a/modules/mol/alg/pymod/scoring.py +++ b/modules/mol/alg/pymod/scoring.py @@ -886,8 +886,11 @@ class Scorer: if trg_ch1 in flat_mapping and trg_ch2 in flat_mapping: mdl_ch1 = flat_mapping[trg_ch1] mdl_ch2 = flat_mapping[trg_ch2] + aln1 = self.mapping.alns[(trg_ch1, mdl_ch1)] + aln2 = self.mapping.alns[(trg_ch2, mdl_ch2)] res = dockq.DockQ(self.model, self.target, mdl_ch1, mdl_ch2, - trg_ch1, trg_ch2) + trg_ch1, trg_ch2, ch1_aln=aln1, + ch2_aln=aln2) if res["nnat"] > 0: self._dockq_interfaces.append((trg_ch1, trg_ch2, mdl_ch1, mdl_ch2)) @@ -897,6 +900,7 @@ class Scorer: # interface which is not covered by mdl... let's run DockQ # with trg as trg/mdl in order to get the native contacts # out + # no need to pass alns as the residue numbers match for sure res = dockq.DockQ(self.target, self.target, trg_ch1, trg_ch2, trg_ch1, trg_ch2) nnat = res["nnat"]