"""
Chain mapping aims to identify a one-to-one relationship between chains in a
reference structure and a model.
"""
import itertools
import copy
import numpy as np
from scipy.special import factorial
from scipy.special import binom # as of Python 3.8, the math module implements
# comb, i.e. n choose k
import ost
from ost import seq
from ost import mol
from ost import geom
from ost.mol.alg import lddt
from ost.mol.alg import qsscore
def _CSel(ent, cnames):
""" Returns view with specified chains
Ensures that quotation marks are around chain names to not confuse
OST query language with weird special characters.
"""
query = "cname=" + ','.join([mol.QueryQuoteName(cname) for cname in cnames])
return ent.Select(query)
class MappingResult:
""" Result object for the chain mapping functions in :class:`ChainMapper`
Constructor is directly called within the functions, no need to construct
such objects yourself.
"""
def __init__(self, target, model, chem_groups, chem_mapping, mapping, alns,
opt_score=None):
self._target = target
self._model = model
self._chem_groups = chem_groups
self._chem_mapping = chem_mapping
self._mapping = mapping
self._alns = alns
self._opt_score = opt_score
@property
def target(self):
""" Target/reference structure, i.e. :attr:`ChainMapper.target`
:type: :class:`ost.mol.EntityView`
"""
return self._target
@property
def model(self):
""" Model structure that gets mapped onto :attr:`~target`
Underwent same processing as :attr:`ChainMapper.target`, i.e.
only contains peptide/nucleotide chains of sufficient size.
:type: :class:`ost.mol.EntityView`
"""
return self._model
@property
def chem_groups(self):
""" Groups of chemically equivalent chains in :attr:`~target`
Same as :attr:`ChainMapper.chem_group`
:class:`list` of :class:`list` of :class:`str` (chain names)
"""
return self._chem_groups
@property
def chem_mapping(self):
""" Assigns chains in :attr:`~model` to :attr:`~chem_groups`.
:class:`list` of :class:`list` of :class:`str` (chain names)
"""
return self._chem_mapping
@property
def mapping(self):
""" Mapping of :attr:`~model` chains onto :attr:`~target`
Exact same shape as :attr:`~chem_groups` but containing the names of the
mapped chains in :attr:`~model`. May contain None for :attr:`~target`
chains that are not covered. No guarantee that all chains in
:attr:`~model` are mapped.
:class:`list` of :class:`list` of :class:`str` (chain names)
"""
return self._mapping
@property
def alns(self):
""" Alignments of mapped chains in :attr:`~target` and :attr:`~model`
Each alignment is accessible with ``alns[(t_chain,m_chain)]``. First
sequence is the sequence of :attr:`target` chain, second sequence the
one from :attr:`~model`. The respective :class:`ost.mol.EntityView` are
attached with :func:`ost.seq.ConstSequenceHandle.AttachView`.
:type: :class:`dict` with key: :class:`tuple` of :class:`str`, value:
:class:`ost.seq.AlignmentHandle`
"""
return self._alns
@property
def opt_score(self):
""" Placeholder property without any guarantee of being set
Different scores get optimized in the various chain mapping algorithms.
Some of them may set their final optimal score in that property.
Consult the documentation of the respective chain mapping algorithm
for more information. Won't be in the return dict of
:func:`JSONSummary`.
"""
return self._opt_score
def GetFlatMapping(self, mdl_as_key=False):
""" Returns flat mapping as :class:`dict` for all mapable chains
:param mdl_as_key: Default is target chain name as key and model chain
name as value. This can be reversed with this flag.
:returns: :class:`dict` with :class:`str` as key/value that describe
one-to-one mapping
"""
flat_mapping = dict()
for trg_chem_group, mdl_chem_group in zip(self.chem_groups,
self.mapping):
for a,b in zip(trg_chem_group, mdl_chem_group):
if a is not None and b is not None:
if mdl_as_key:
flat_mapping[b] = a
else:
flat_mapping[a] = b
return flat_mapping
def JSONSummary(self):
""" Returns JSON serializable summary of results
"""
json_dict = dict()
json_dict["chem_groups"] = self.chem_groups
json_dict["mapping"] = self.mapping
json_dict["flat_mapping"] = self.GetFlatMapping()
json_dict["alns"] = list()
for aln in self.alns.values():
trg_seq = aln.GetSequence(0)
mdl_seq = aln.GetSequence(1)
aln_dict = {"trg_ch": trg_seq.GetName(), "trg_seq": str(trg_seq),
"mdl_ch": mdl_seq.GetName(), "mdl_seq": str(mdl_seq)}
json_dict["alns"].append(aln_dict)
return json_dict
class ReprResult:
""" Result object for :func:`ChainMapper.GetRepr`
Constructor is directly called within the function, no need to construct
such objects yourself.
:param lDDT: lDDT for this mapping. Depends on how you call
:func:`ChainMapper.GetRepr` whether this is backbone only or
full atom lDDT.
:type lDDT: :class:`float`
:param substructure: The full substructure for which we searched for a
representation
:type substructure: :class:`ost.mol.EntityView`
:param ref_view: View pointing to the same underlying entity as
*substructure* but only contains the stuff that is mapped
:type ref_view: :class:`mol.EntityView`
:param mdl_view: The matching counterpart in model
:type mdl_view: :class:`mol.EntityView`
"""
def __init__(self, lDDT, substructure, ref_view, mdl_view):
self._lDDT = lDDT
self._substructure = substructure
assert(len(ref_view.residues) == len(mdl_view.residues))
self._ref_view = ref_view
self._mdl_view = mdl_view
# lazily evaluated attributes
self._ref_bb_pos = None
self._mdl_bb_pos = None
self._ref_full_bb_pos = None
self._mdl_full_bb_pos = None
self._transform = None
self._superposed_mdl_bb_pos = None
self._bb_rmsd = None
self._gdt_8 = None
self._gdt_4 = None
self._gdt_2 = None
self._gdt_1 = None
self._ost_query = None
self._flat_mapping = None
self._inconsistent_residues = None
@property
def lDDT(self):
""" lDDT of representation result
Depends on how you call :func:`ChainMapper.GetRepr` whether this is
backbone only or full atom lDDT.
:type: :class:`float`
"""
return self._lDDT
@property
def substructure(self):
""" The full substructure for which we searched for a
representation
:type: :class:`ost.mol.EntityView`
"""
return self._substructure
@property
def ref_view(self):
""" View which contains the mapped subset of :attr:`substructure`
:type: :class:`ost.mol.EntityView`
"""
return self._ref_view
@property
def mdl_view(self):
""" The :attr:`ref_view` representation in the model
:type: :class:`ost.mol.EntityView`
"""
return self._mdl_view
@property
def ref_residues(self):
""" The reference residues
:type: class:`mol.ResidueViewList`
"""
return self.ref_view.residues
@property
def mdl_residues(self):
""" The model residues
:type: :class:`mol.ResidueViewList`
"""
return self.mdl_view.residues
@property
def inconsistent_residues(self):
""" A list of mapped residue whose names do not match (eg. ALA in the
reference and LEU in the model).
The mismatches are reported as a tuple of :class:`~ost.mol.ResidueView`
(reference, model), or as an empty list if all the residue names match.
:type: :class:`list`
"""
if self._inconsistent_residues is None:
self._inconsistent_residues = self._GetInconsistentResidues(
self.ref_residues, self.mdl_residues)
return self._inconsistent_residues
@property
def ref_bb_pos(self):
""" Representative backbone positions for reference residues.
Thats CA positions for peptides and C3' positions for Nucleotides.
:type: :class:`geom.Vec3List`
"""
if self._ref_bb_pos is None:
self._ref_bb_pos = self._GetBBPos(self.ref_residues)
return self._ref_bb_pos
@property
def mdl_bb_pos(self):
""" Representative backbone positions for model residues.
Thats CA positions for peptides and C3' positions for Nucleotides.
:type: :class:`geom.Vec3List`
"""
if self._mdl_bb_pos is None:
self._mdl_bb_pos = self._GetBBPos(self.mdl_residues)
return self._mdl_bb_pos
@property
def ref_full_bb_pos(self):
""" Representative backbone positions for reference residues.
Thats N, CA and C positions for peptides and O5', C5', C4', C3', O3'
positions for Nucleotides.
:type: :class:`geom.Vec3List`
"""
if self._ref_full_bb_pos is None:
self._ref_full_bb_pos = self._GetFullBBPos(self.ref_residues)
return self._ref_full_bb_pos
@property
def mdl_full_bb_pos(self):
""" Representative backbone positions for reference residues.
Thats N, CA and C positions for peptides and O5', C5', C4', C3', O3'
positions for Nucleotides.
:type: :class:`geom.Vec3List`
"""
if self._mdl_full_bb_pos is None:
self._mdl_full_bb_pos = self._GetFullBBPos(self.mdl_residues)
return self._mdl_full_bb_pos
@property
def transform(self):
""" Transformation to superpose mdl residues onto ref residues
Superposition computed as minimal RMSD superposition on
:attr:`ref_bb_pos` and :attr:`mdl_bb_pos`. If number of positions is
smaller 3, the full_bb_pos equivalents are used instead.
:type: :class:`ost.geom.Mat4`
"""
if self._transform is None:
if len(self.mdl_bb_pos) < 3:
self._transform = _GetTransform(self.mdl_full_bb_pos,
self.ref_full_bb_pos, False)
else:
self._transform = _GetTransform(self.mdl_bb_pos,
self.ref_bb_pos, False)
return self._transform
@property
def superposed_mdl_bb_pos(self):
""" :attr:`mdl_bb_pos` with :attr:`transform applied`
:type: :class:`geom.Vec3List`
"""
if self._superposed_mdl_bb_pos is None:
self._superposed_mdl_bb_pos = geom.Vec3List(self.mdl_bb_pos)
self._superposed_mdl_bb_pos.ApplyTransform(self.transform)
return self._superposed_mdl_bb_pos
@property
def bb_rmsd(self):
""" RMSD between :attr:`ref_bb_pos` and :attr:`superposed_mdl_bb_pos`
:type: :class:`float`
"""
if self._bb_rmsd is None:
self._bb_rmsd = self.ref_bb_pos.GetRMSD(self.superposed_mdl_bb_pos)
return self._bb_rmsd
@property
def gdt_8(self):
""" GDT with one single threshold: 8.0
:type: :class:`float`
"""
if self._gdt_8 is None:
self._gdt_8 = self.ref_bb_pos.GetGDT(self.superposed_mdl_bb_pos, 8.0)
return self._gdt_8
@property
def gdt_4(self):
""" GDT with one single threshold: 4.0
:type: :class:`float`
"""
if self._gdt_4 is None:
self._gdt_4 = self.ref_bb_pos.GetGDT(self.superposed_mdl_bb_pos, 4.0)
return self._gdt_4
@property
def gdt_2(self):
""" GDT with one single threshold: 2.0
:type: :class:`float`
"""
if self._gdt_2 is None:
self._gdt_2 = self.ref_bb_pos.GetGDT(self.superposed_mdl_bb_pos, 2.0)
return self._gdt_2
@property
def gdt_1(self):
""" GDT with one single threshold: 1.0
:type: :class:`float`
"""
if self._gdt_1 is None:
self._gdt_1 = self.ref_bb_pos.GetGDT(self.superposed_mdl_bb_pos, 1.0)
return self._gdt_1
@property
def ost_query(self):
""" query for mdl residues in OpenStructure query language
Repr can be selected as ``full_mdl.Select(ost_query)``
Returns invalid query if residue numbers have insertion codes.
:type: :class:`str`
"""
if self._ost_query is None:
chain_rnums = dict()
for r in self.mdl_residues:
chname = r.GetChain().GetName()
rnum = r.GetNumber().GetNum()
if chname not in chain_rnums:
chain_rnums[chname] = list()
chain_rnums[chname].append(str(rnum))
chain_queries = list()
for k,v in chain_rnums.items():
q = f"(cname={mol.QueryQuoteName(k)} and "
q += f"rnum={','.join(v)})"
chain_queries.append(q)
self._ost_query = " or ".join(chain_queries)
return self._ost_query
def JSONSummary(self):
""" Returns JSON serializable summary of results
"""
json_dict = dict()
json_dict["lDDT"] = self.lDDT
json_dict["ref_residues"] = [r.GetQualifiedName() for r in \
self.ref_residues]
json_dict["mdl_residues"] = [r.GetQualifiedName() for r in \
self.mdl_residues]
json_dict["transform"] = list(self.transform.data)
json_dict["bb_rmsd"] = self.bb_rmsd
json_dict["gdt_8"] = self.gdt_8
json_dict["gdt_4"] = self.gdt_4
json_dict["gdt_2"] = self.gdt_2
json_dict["gdt_1"] = self.gdt_1
json_dict["ost_query"] = self.ost_query
json_dict["flat_mapping"] = self.GetFlatChainMapping()
return json_dict
def GetFlatChainMapping(self, mdl_as_key=False):
""" Returns flat mapping of all chains in the representation
:param mdl_as_key: Default is target chain name as key and model chain
name as value. This can be reversed with this flag.
:returns: :class:`dict` with :class:`str` as key/value that describe
one-to-one mapping
"""
flat_mapping = dict()
for trg_res, mdl_res in zip(self.ref_residues, self.mdl_residues):
if mdl_as_key:
flat_mapping[mdl_res.chain.name] = trg_res.chain.name
else:
flat_mapping[trg_res.chain.name] = mdl_res.chain.name
return flat_mapping
def _GetFullBBPos(self, residues):
""" Helper to extract full backbone positions
"""
exp_pep_atoms = ["N", "CA", "C"]
exp_nuc_atoms = ["\"O5'\"", "\"C5'\"", "\"C4'\"", "\"C3'\"", "\"O3'\""]
bb_pos = geom.Vec3List()
for r in residues:
if r.GetChemType() == mol.ChemType.NUCLEOTIDES:
exp_atoms = exp_nuc_atoms
elif r.GetChemType() == mol.ChemType.AMINOACIDS:
exp_atoms = exp_pep_atoms
else:
raise RuntimeError("Something terrible happened... RUN...")
for aname in exp_atoms:
a = r.FindAtom(aname)
if not a.IsValid():
raise RuntimeError("Something terrible happened... "
"RUN...")
bb_pos.append(a.GetPos())
return bb_pos
def _GetBBPos(self, residues):
""" Helper to extract single representative position for each residue
"""
bb_pos = geom.Vec3List()
for r in residues:
at = r.FindAtom("CA")
if not at.IsValid():
at = r.FindAtom("C3'")
if not at.IsValid():
raise RuntimeError("Something terrible happened... RUN...")
bb_pos.append(at.GetPos())
return bb_pos
def _GetInconsistentResidues(self, ref_residues, mdl_residues):
""" Helper to extract a list of inconsistent residues.
"""
if len(ref_residues) != len(mdl_residues):
raise ValueError("Something terrible happened... Reference and "
"model lengths differ... RUN...")
inconsistent_residues = list()
for ref_residue, mdl_residue in zip(ref_residues, mdl_residues):
if ref_residue.name != mdl_residue.name:
inconsistent_residues.append((ref_residue, mdl_residue))
return inconsistent_residues
class ChainMapper:
""" Class to compute chain mappings
All algorithms are performed on processed structures which fulfill
criteria as given in constructor arguments (*min_pep_length*,
"min_nuc_length") and only contain residues which have all required backbone
atoms. for peptide residues thats N, CA, C and CB (no CB for GLY), for
nucleotide residues thats O5', C5', C4', C3' and O3'.
Chain mapping is a three step process:
* Group chemically identical chains in *target* using pairwise
alignments that are either computed with Needleman-Wunsch (NW) or
simply derived from residue numbers (*resnum_alignments* flag).
In case of NW, *pep_subst_mat*, *pep_gap_open* and *pep_gap_ext*
and their nucleotide equivalents are relevant. Two chains are
considered identical if they fulfill the thresholds given by
*pep_seqid_thr*, *pep_gap_thr*, their nucleotide equivalents
respectively. The grouping information is available as
attributes of this class.
* Map chains in an input model to these groups. Generating alignments
and the similarity criteria are the same as above. You can either
get the group mapping with :func:`GetChemMapping` or directly call
one of the full fletched one-to-one chain mapping functions which
execute that step internally.
* Obtain one-to-one mapping for chains in an input model and
*target* with one of the available mapping functions. Just to get an
idea of complexity. If *target* and *model* are octamers, there are
``8! = 40320`` possible chain mappings.
:param target: Target structure onto which models are mapped.
Computations happen on a selection only containing
polypeptides and polynucleotides.
:type target: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
:param resnum_alignments: Use residue numbers instead of
Needleman-Wunsch to compute pairwise
alignments. Relevant for :attr:`~chem_groups`
and related attributes.
:type resnum_alignments: :class:`bool`
:param pep_seqid_thr: Threshold used to decide when two chains are
identical. 95 percent tolerates the few mutations
crystallographers like to do.
:type pep_seqid_thr: :class:`float`
:param pep_gap_thr: Additional threshold to avoid gappy alignments with
high seqid. By default this is disabled (set to 1.0).
This threshold checks for a maximum allowed fraction
of gaps in any of the two sequences after stripping
terminal gaps. The reason for not just normalizing
seqid by the longer sequence is that one sequence
might be a perfect subsequence of the other but only
cover half of it.
:type pep_gap_thr: :class:`float`
:param nuc_seqid_thr: Nucleotide equivalent for *pep_seqid_thr*
:type nuc_seqid_thr: :class:`float`
:param nuc_gap_thr: Nucleotide equivalent for *nuc_gap_thr*
:type nuc_gap_thr: :class:`float`
:param pep_subst_mat: Substitution matrix to align peptide sequences,
irrelevant if *resnum_alignments* is True,
defaults to seq.alg.BLOSUM62
:type pep_subst_mat: :class:`ost.seq.alg.SubstWeightMatrix`
:param pep_gap_open: Gap open penalty to align peptide sequences,
irrelevant if *resnum_alignments* is True
:type pep_gap_open: :class:`int`
:param pep_gap_ext: Gap extension penalty to align peptide sequences,
irrelevant if *resnum_alignments* is True
:type pep_gap_ext: :class:`int`
:param nuc_subst_mat: Nucleotide equivalent for *pep_subst_mat*,
defaults to seq.alg.NUC44
:type nuc_subst_mat: :class:`ost.seq.alg.SubstWeightMatrix`
:param nuc_gap_open: Nucleotide equivalent for *pep_gap_open*
:type nuc_gap_open: :class:`int`
:param nuc_gap_ext: Nucleotide equivalent for *pep_gap_ext*
:type nuc_gap_ext: :class:`int`
:param min_pep_length: Minimal number of residues for a peptide chain to be
considered in target and in models.
:type min_pep_length: :class:`int`
:param min_nuc_length: Minimal number of residues for a nucleotide chain to be
considered in target and in models.
:type min_nuc_length: :class:`int`
:param n_max_naive: Max possible chain mappings that are enumerated in
:func:`~GetNaivelDDTMapping` /
:func:`~GetDecomposerlDDTMapping`. A
:class:`RuntimeError` is raised in case of bigger
complexity.
:type n_max_naive: :class:`int`
"""
def __init__(self, target, resnum_alignments=False,
pep_seqid_thr = 95., pep_gap_thr = 1.0,
nuc_seqid_thr = 95., nuc_gap_thr = 1.0,
pep_subst_mat = seq.alg.BLOSUM62, pep_gap_open = -11,
pep_gap_ext = -1, nuc_subst_mat = seq.alg.NUC44,
nuc_gap_open = -4, nuc_gap_ext = -4,
min_pep_length = 10, min_nuc_length = 4,
n_max_naive = 1e8):
# attributes
self.resnum_alignments = resnum_alignments
self.pep_seqid_thr = pep_seqid_thr
self.pep_gap_thr = pep_gap_thr
self.nuc_seqid_thr = nuc_seqid_thr
self.nuc_gap_thr = nuc_gap_thr
self.min_pep_length = min_pep_length
self.min_nuc_length = min_nuc_length
self.n_max_naive = n_max_naive
# lazy computed attributes
self._chem_groups = None
self._chem_group_alignments = None
self._chem_group_ref_seqs = None
self._chem_group_types = None
# helper class to generate pairwise alignments
self.aligner = _Aligner(resnum_aln = resnum_alignments,
pep_subst_mat = pep_subst_mat,
pep_gap_open = pep_gap_open,
pep_gap_ext = pep_gap_ext,
nuc_subst_mat = nuc_subst_mat,
nuc_gap_open = nuc_gap_open,
nuc_gap_ext = nuc_gap_ext)
# target structure preprocessing
self._target, self._polypep_seqs, self._polynuc_seqs = \
self.ProcessStructure(target)
@property
def target(self):
"""Target structure that only contains peptides/nucleotides
Contains only residues that have the backbone representatives
(CA for peptide and C3' for nucleotides) to avoid ATOMSEQ alignment
inconsistencies when switching between all atom and backbone only
representations.
:type: :class:`ost.mol.EntityView`
"""
return self._target
@property
def polypep_seqs(self):
"""Sequences of peptide chains in :attr:`~target`
Respective :class:`EntityView` from *target* for each sequence s are
available as ``s.GetAttachedView()``
:type: :class:`ost.seq.SequenceList`
"""
return self._polypep_seqs
@property
def polynuc_seqs(self):
"""Sequences of nucleotide chains in :attr:`~target`
Respective :class:`EntityView` from *target* for each sequence s are
available as ``s.GetAttachedView()``
:type: :class:`ost.seq.SequenceList`
"""
return self._polynuc_seqs
@property
def chem_groups(self):
"""Groups of chemically equivalent chains in :attr:`~target`
First chain in group is the one with longest sequence.
:getter: Computed on first use (cached)
:type: :class:`list` of :class:`list` of :class:`str` (chain names)
"""
if self._chem_groups is None:
self._chem_groups = list()
for a in self.chem_group_alignments:
self._chem_groups.append([s.GetName() for s in a.sequences])
return self._chem_groups
@property
def chem_group_alignments(self):
"""MSA for each group in :attr:`~chem_groups`
Sequences in MSAs exhibit same order as in :attr:`~chem_groups` and
have the respective :class:`ost.mol.EntityView` from *target* attached.
:getter: Computed on first use (cached)
:type: :class:`ost.seq.AlignmentList`
"""
if self._chem_group_alignments is None:
self._chem_group_alignments, self._chem_group_types = \
_GetChemGroupAlignments(self.polypep_seqs, self.polynuc_seqs,
self.aligner,
pep_seqid_thr=self.pep_seqid_thr,
pep_gap_thr=self.pep_gap_thr,
nuc_seqid_thr=self.nuc_seqid_thr,
nuc_gap_thr=self.nuc_gap_thr)
return self._chem_group_alignments
@property
def chem_group_ref_seqs(self):
"""Reference (longest) sequence for each group in :attr:`~chem_groups`
Respective :class:`EntityView` from *target* for each sequence s are
available as ``s.GetAttachedView()``
:getter: Computed on first use (cached)
:type: :class:`ost.seq.SequenceList`
"""
if self._chem_group_ref_seqs is None:
self._chem_group_ref_seqs = seq.CreateSequenceList()
for a in self.chem_group_alignments:
s = seq.CreateSequence(a.GetSequence(0).GetName(),
a.GetSequence(0).GetGaplessString())
s.AttachView(a.GetSequence(0).GetAttachedView())
self._chem_group_ref_seqs.AddSequence(s)
return self._chem_group_ref_seqs
@property
def chem_group_types(self):
"""ChemType of each group in :attr:`~chem_groups`
Specifying if groups are poly-peptides/nucleotides, i.e.
:class:`ost.mol.ChemType.AMINOACIDS` or
:class:`ost.mol.ChemType.NUCLEOTIDES`
:getter: Computed on first use (cached)
:type: :class:`list` of :class:`ost.mol.ChemType`
"""
if self._chem_group_types is None:
self._chem_group_alignments, self._chem_group_types = \
_GetChemGroupAlignments(self.polypep_seqs, self.polynuc_seqs,
self.aligner,
pep_seqid_thr=self.pep_seqid_thr,
pep_gap_thr=self.pep_gap_thr,
nuc_seqid_thr=self.nuc_seqid_thr,
nuc_gap_thr=self.nuc_gap_thr)
return self._chem_group_types
def GetChemMapping(self, model):
"""Maps sequences in *model* to chem_groups of target
:param model: Model from which to extract sequences, a
selection that only includes peptides and nucleotides
is performed and returned along other results.
:type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
:returns: Tuple with two lists of length `len(self.chem_groups)` and
an :class:`ost.mol.EntityView` representing *model*:
1) Each element is a :class:`list` with mdl chain names that
map to the chem group at that position.
2) Each element is a :class:`ost.seq.AlignmentList` aligning
these mdl chain sequences to the chem group ref sequences.
3) A selection of *model* that only contains polypeptides and
polynucleotides whose ATOMSEQ exactly matches the sequence
info in the returned alignments.
"""
mdl, mdl_pep_seqs, mdl_nuc_seqs = self.ProcessStructure(model)
mapping = [list() for x in self.chem_groups]
alns = [seq.AlignmentList() for x in self.chem_groups]
for s in mdl_pep_seqs:
idx, aln = _MapSequence(self.chem_group_ref_seqs,
self.chem_group_types,
s, mol.ChemType.AMINOACIDS,
self.aligner)
if idx is not None:
mapping[idx].append(s.GetName())
alns[idx].append(aln)
for s in mdl_nuc_seqs:
idx, aln = _MapSequence(self.chem_group_ref_seqs,
self.chem_group_types,
s, mol.ChemType.NUCLEOTIDES,
self.aligner)
if idx is not None:
mapping[idx].append(s.GetName())
alns[idx].append(aln)
return (mapping, alns, mdl)
def GetlDDTMapping(self, model, inclusion_radius=15.0,
thresholds=[0.5, 1.0, 2.0, 4.0], strategy="naive",
steep_opt_rate = None, block_seed_size = 5,
block_blocks_per_chem_group = 5,
chem_mapping_result = None):
""" Identify chain mapping by optimizing lDDT score
Maps *model* chain sequences to :attr:`~chem_groups` and find mapping
based on backbone only lDDT score (CA for amino acids C3' for
Nucleotides).
Either performs a naive search, i.e. enumerate all possible mappings or
executes a greedy strategy that tries to identify a (close to) optimal
mapping in an iterative way by starting from a start mapping (seed). In
each iteration, the one-to-one mapping that leads to highest increase
in number of conserved contacts is added with the additional requirement
that this added mapping must have non-zero interface counts towards the
already mapped chains. So basically we're "growing" the mapped structure
by only adding connected stuff.
The available strategies:
* **naive**: Enumerates all possible mappings and returns best
* **greedy_fast**: perform all vs. all single chain lDDTs within the
respective ref/mdl chem groups. The mapping with highest number of
conserved contacts is selected as seed for greedy extension
* **greedy_full**: try multiple seeds for greedy extension, i.e. try
all ref/mdl chain combinations within the respective chem groups and
retain the mapping leading to the best lDDT.
* **greedy_block**: try multiple seeds for greedy extension, i.e. try
all ref/mdl chain combinations within the respective chem groups and
extend them to *block_seed_size*. *block_blocks_per_chem_group*
for each chem group are selected for exhaustive extension.
Sets :attr:`MappingResult.opt_score` in case of no trivial one-to-one
mapping.
:param model: Model to map
:type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
:param inclusion_radius: Inclusion radius for lDDT
:type inclusion_radius: :class:`float`
:param thresholds: Thresholds for lDDT
:type thresholds: :class:`list` of :class:`float`
:param strategy: Strategy to find mapping. Must be in ["naive",
"greedy_fast", "greedy_full", "greedy_block"]
:type strategy: :class:`str`
:param steep_opt_rate: Only relevant for greedy strategies.
If set, every *steep_opt_rate* mappings, a simple
optimization is executed with the goal of
avoiding local minima. The optimization
iteratively checks all possible swaps of mappings
within their respective chem groups and accepts
swaps that improve lDDT score. Iteration stops as
soon as no improvement can be achieved anymore.
:type steep_opt_rate: :class:`int`
:param block_seed_size: Param for *greedy_block* strategy - Initial seeds
are extended by that number of chains.
:type block_seed_size: :class:`int`
:param block_blocks_per_chem_group: Param for *greedy_block* strategy -
Number of blocks per chem group that
are extended in an initial search
for high scoring local solutions.
:type block_blocks_per_chem_group: :class:`int`
:param chem_mapping_result: Pro param. The result of
:func:`~GetChemMapping` where you provided
*model*. If set, *model* parameter is not
used.
:type chem_mapping_result: :class:`tuple`
:returns: A :class:`MappingResult`
"""
strategies = ["naive", "greedy_fast", "greedy_full", "greedy_block"]
if strategy not in strategies:
raise RuntimeError(f"Strategy must be in {strategies}")
if chem_mapping_result is None:
chem_mapping, chem_group_alns, mdl = self.GetChemMapping(model)
else:
chem_mapping, chem_group_alns, mdl = chem_mapping_result
ref_mdl_alns = _GetRefMdlAlns(self.chem_groups,
self.chem_group_alignments,
chem_mapping,
chem_group_alns)
# check for the simplest case
one_to_one = _CheckOneToOneMapping(self.chem_groups, chem_mapping)
if one_to_one is not None:
alns = dict()
for ref_group, mdl_group in zip(self.chem_groups, one_to_one):
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
aln.AttachView(0, _CSel(self.target, [ref_ch]))
aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
one_to_one, alns)
mapping = None
opt_lddt = None
if strategy == "naive":
mapping, opt_lddt = _lDDTNaive(self.target, mdl, inclusion_radius,
thresholds, self.chem_groups,
chem_mapping, ref_mdl_alns,
self.n_max_naive)
else:
# its one of the greedy strategies - setup greedy searcher
the_greed = _lDDTGreedySearcher(self.target, mdl, self.chem_groups,
chem_mapping, ref_mdl_alns,
inclusion_radius=inclusion_radius,
thresholds=thresholds,
steep_opt_rate=steep_opt_rate)
if strategy == "greedy_fast":
mapping = _lDDTGreedyFast(the_greed)
elif strategy == "greedy_full":
mapping = _lDDTGreedyFull(the_greed)
elif strategy == "greedy_block":
mapping = _lDDTGreedyBlock(the_greed, block_seed_size,
block_blocks_per_chem_group)
# cached => lDDT computation is fast here
opt_lddt = the_greed.lDDT(self.chem_groups, mapping)
alns = dict()
for ref_group, mdl_group in zip(self.chem_groups, mapping):
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
aln.AttachView(0, _CSel(self.target, [ref_ch]))
aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
mapping, alns, opt_score = opt_lddt)
def GetQSScoreMapping(self, model, contact_d = 12.0, strategy = "naive",
block_seed_size = 5, block_blocks_per_chem_group = 5,
steep_opt_rate = None, chem_mapping_result = None,
greedy_prune_contact_map = False):
""" Identify chain mapping based on QSScore
Scoring is based on CA/C3' positions which are present in all chains of
a :attr:`chem_groups` as well as the *model* chains which are mapped to
that respective chem group.
The following strategies are available:
* **naive**: Naively iterate all possible mappings and return best based
on QS score.
* **greedy_fast**: perform all vs. all single chain lDDTs within the
respective ref/mdl chem groups. The mapping with highest number of
conserved contacts is selected as seed for greedy extension.
Extension is based on QS-score.
* **greedy_full**: try multiple seeds for greedy extension, i.e. try
all ref/mdl chain combinations within the respective chem groups and
retain the mapping leading to the best QS-score.
* **greedy_block**: try multiple seeds for greedy extension, i.e. try
all ref/mdl chain combinations within the respective chem groups and
extend them to *block_seed_size*. *block_blocks_per_chem_group*
for each chem group are selected for exhaustive extension.
Sets :attr:`MappingResult.opt_score` in case of no trivial one-to-one
mapping.
:param model: Model to map
:type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
:param contact_d: Max distance between two residues to be considered as
contact in qs scoring
:type contact_d: :class:`float`
:param strategy: Strategy for sampling, must be in ["naive"]
:type strategy: :class:`str`
:param chem_mapping_result: Pro param. The result of
:func:`~GetChemMapping` where you provided
*model*. If set, *model* parameter is not
used.
:type chem_mapping_result: :class:`tuple`
:param greedy_prune_contact_map: Relevant for all strategies that use
greedy extensions. If True, only chains
with at least 3 contacts (8A CB
distance) towards already mapped chains
in trg/mdl are considered for
extension. All chains that give a
potential non-zero QS-score increase
are used otherwise (at least one
contact within 12A). The consequence
is reduced runtime and usually no
real reduction in accuracy.
:returns: A :class:`MappingResult`
"""
strategies = ["naive", "greedy_fast", "greedy_full", "greedy_block"]
if strategy not in strategies:
raise RuntimeError(f"strategy must be {strategies}")
if chem_mapping_result is None:
chem_mapping, chem_group_alns, mdl = self.GetChemMapping(model)
else:
chem_mapping, chem_group_alns, mdl = chem_mapping_result
ref_mdl_alns = _GetRefMdlAlns(self.chem_groups,
self.chem_group_alignments,
chem_mapping,
chem_group_alns)
# check for the simplest case
one_to_one = _CheckOneToOneMapping(self.chem_groups, chem_mapping)
if one_to_one is not None:
alns = dict()
for ref_group, mdl_group in zip(self.chem_groups, one_to_one):
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
aln.AttachView(0, _CSel(self.target, [ref_ch]))
aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
one_to_one, alns)
mapping = None
opt_qsscore = None
if strategy == "naive":
mapping, opt_qsscore = _QSScoreNaive(self.target, mdl,
self.chem_groups,
chem_mapping, ref_mdl_alns,
contact_d, self.n_max_naive)
else:
# its one of the greedy strategies - setup greedy searcher
the_greed = _QSScoreGreedySearcher(self.target, mdl,
self.chem_groups,
chem_mapping, ref_mdl_alns,
contact_d = contact_d,
steep_opt_rate=steep_opt_rate,
greedy_prune_contact_map = greedy_prune_contact_map)
if strategy == "greedy_fast":
mapping = _QSScoreGreedyFast(the_greed)
elif strategy == "greedy_full":
mapping = _QSScoreGreedyFull(the_greed)
elif strategy == "greedy_block":
mapping = _QSScoreGreedyBlock(the_greed, block_seed_size,
block_blocks_per_chem_group)
# cached => QSScore computation is fast here
opt_qsscore = the_greed.Score(mapping, check=False)
alns = dict()
for ref_group, mdl_group in zip(self.chem_groups, mapping):
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
aln.AttachView(0, _CSel(self.target, [ref_ch]))
aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
mapping, alns, opt_score = opt_qsscore)
def GetRigidMapping(self, model, strategy = "greedy_single_gdtts",
single_chain_gdtts_thresh=0.4, subsampling=None,
first_complete=False, iterative_superposition=False,
chem_mapping_result = None):
"""Identify chain mapping based on rigid superposition
Superposition and scoring is based on CA/C3' positions which are present
in all chains of a :attr:`chem_groups` as well as the *model*
chains which are mapped to that respective chem group.
Transformations to superpose *model* onto :attr:`ChainMapper.target`
are estimated using all possible combinations of target and model chains
within the same chem groups and build the basis for further extension.
There are four extension strategies:
* **greedy_single_gdtts**: Iteratively add the model/target chain pair
that adds the most conserved contacts based on the GDT-TS metric
(Number of CA/C3' atoms within [8, 4, 2, 1] Angstrom). The mapping
with highest GDT-TS score is returned. However, that mapping is not
guaranteed to be complete (see *single_chain_gdtts_thresh*).
* **greedy_iterative_gdtts**: Same as greedy_single_gdtts except that
the transformation gets updated with each added chain pair.
* **greedy_single_rmsd**: Conceptually similar to greedy_single_gdtts
but the added chain pairs are the ones with lowest RMSD.
The mapping with lowest overall RMSD gets returned.
*single_chain_gdtts_thresh* is only applied to derive the initial
transformations. After that, the minimal RMSD chain pair gets
iteratively added without applying any threshold.
* **greedy_iterative_rmsd**: Same as greedy_single_rmsd exept that
the transformation gets updated with each added chain pair.
*single_chain_gdtts_thresh* is only applied to derive the initial
transformations. After that, the minimal RMSD chain pair gets
iteratively added without applying any threshold.
:param model: Model to map
:type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
:param strategy: Strategy to extend mappings from initial transforms,
see description above. Must be in ["greedy_single",
"greedy_iterative", "greedy_iterative_rmsd"]
:type strategy: :class:`str`
:param single_chain_gdtts_thresh: Minimal GDT-TS score for model/target
chain pair to be added to mapping.
Mapping extension for a given
transform stops when no pair fulfills
this threshold, potentially leading to
an incomplete mapping.
:type single_chain_gdtts_thresh: :class:`float`
:param subsampling: If given, only use an equally distributed subset
of all CA/C3' positions for superposition/scoring.
:type subsampling: :class:`int`
:param first_complete: Avoid full enumeration and return first found
mapping that covers all model chains or all
target chains. Has no effect on
greedy_iterative_rmsd strategy.
:type first_complete: :class:`bool`
:param iterative_superposition: Whether to compute inital
transformations with
:func:`ost.mol.alg.IterativeSuperposeSVD`
as oposed to
:func:`ost.mol.alg.SuperposeSVD`
:type iterative_superposition: :class:`bool`
:param chem_mapping_result: Pro param. The result of
:func:`~GetChemMapping` where you provided
*model*. If set, *model* parameter is not
used.
:type chem_mapping_result: :class:`tuple`
:returns: A :class:`MappingResult`
"""
strategies = ["greedy_single_gdtts", "greedy_iterative_gdtts",
"greedy_single_rmsd", "greedy_iterative_rmsd"]
if strategy not in strategies:
raise RuntimeError(f"strategy must be {strategies}")
if chem_mapping_result is None:
chem_mapping, chem_group_alns, mdl = self.GetChemMapping(model)
else:
chem_mapping, chem_group_alns, mdl = chem_mapping_result
ref_mdl_alns = _GetRefMdlAlns(self.chem_groups,
self.chem_group_alignments,
chem_mapping,
chem_group_alns)
# check for the simplest case
one_to_one = _CheckOneToOneMapping(self.chem_groups, chem_mapping)
if one_to_one is not None:
alns = dict()
for ref_group, mdl_group in zip(self.chem_groups, one_to_one):
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
aln.AttachView(0, _CSel(self.target, [ref_ch]))
aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
one_to_one, alns)
trg_group_pos, mdl_group_pos = _GetRefPos(self.target, mdl,
self.chem_group_alignments,
chem_group_alns,
max_pos = subsampling)
# get transforms of any mdl chain onto any trg chain in same chem group
# that fulfills gdtts threshold
initial_transforms = list()
initial_mappings = list()
for trg_pos, trg_chains, mdl_pos, mdl_chains in zip(trg_group_pos,
self.chem_groups,
mdl_group_pos,
chem_mapping):
for t_pos, t in zip(trg_pos, trg_chains):
for m_pos, m in zip(mdl_pos, mdl_chains):
if len(t_pos) >= 3 and len(m_pos) >= 3:
transform = _GetTransform(m_pos, t_pos,
iterative_superposition)
t_m_pos = geom.Vec3List(m_pos)
t_m_pos.ApplyTransform(transform)
gdt = t_pos.GetGDTTS(t_m_pos)
if gdt >= single_chain_gdtts_thresh:
initial_transforms.append(transform)
initial_mappings.append((t,m))
if strategy == "greedy_single_gdtts":
mapping = _SingleRigidGDTTS(initial_transforms, initial_mappings,
self.chem_groups, chem_mapping,
trg_group_pos, mdl_group_pos,
single_chain_gdtts_thresh,
iterative_superposition, first_complete,
len(self.target.chains),
len(mdl.chains))
elif strategy == "greedy_iterative_gdtts":
mapping = _IterativeRigidGDTTS(initial_transforms, initial_mappings,
self.chem_groups, chem_mapping,
trg_group_pos, mdl_group_pos,
single_chain_gdtts_thresh,
iterative_superposition,
first_complete,
len(self.target.chains),
len(mdl.chains))
elif strategy == "greedy_single_rmsd":
mapping = _SingleRigidRMSD(initial_transforms, initial_mappings,
self.chem_groups, chem_mapping,
trg_group_pos, mdl_group_pos,
iterative_superposition)
elif strategy == "greedy_iterative_rmsd":
mapping = _IterativeRigidRMSD(initial_transforms, initial_mappings,
self.chem_groups, chem_mapping,
trg_group_pos, mdl_group_pos,
iterative_superposition)
# translate mapping format and return
final_mapping = list()
for ref_chains in self.chem_groups:
mapped_mdl_chains = list()
for ref_ch in ref_chains:
if ref_ch in mapping:
mapped_mdl_chains.append(mapping[ref_ch])
else:
mapped_mdl_chains.append(None)
final_mapping.append(mapped_mdl_chains)
alns = dict()
for ref_group, mdl_group in zip(self.chem_groups, final_mapping):
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
aln.AttachView(0, _CSel(self.target, [ref_ch]))
aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
final_mapping, alns)
def GetMapping(self, model, n_max_naive = 40320):
""" Convenience function to get mapping with currently preferred method
If number of possible chain mappings is <= *n_max_naive*, a naive
QS-score mapping is performed and optimal QS-score is guaranteed.
For anything else, a QS-score mapping with the greedy_full strategy is
performed (greedy_prune_contact_map = True). The default for
*n_max_naive* of 40320 corresponds to an octamer (8!=40320). A
structure with stoichiometry A6B2 would be 6!*2!=1440 etc.
"""
chem_mapping_res = self.GetChemMapping(model)
if _NMappingsWithin(self.chem_groups, chem_mapping_res[0], n_max_naive):
return self.GetQSScoreMapping(model, strategy="naive",
chem_mapping_result=chem_mapping_res)
else:
return self.GetQSScoreMapping(model, strategy="greedy_full",
greedy_prune_contact_map=True,
chem_mapping_result=chem_mapping_res)
def GetRepr(self, substructure, model, topn=1, inclusion_radius=15.0,
thresholds=[0.5, 1.0, 2.0, 4.0], bb_only=False,
only_interchain=False, chem_mapping_result = None,
global_mapping = None):
""" Identify *topn* representations of *substructure* in *model*
*substructure* defines a subset of :attr:`~target` for which one
wants the *topn* representations in *model*. Representations are scored
and sorted by lDDT.
:param substructure: A :class:`ost.mol.EntityView` which is a subset of
:attr:`~target`. Should be selected with the
OpenStructure query language. Example: if you're
interested in residues with number 42,43 and 85 in
chain A:
``substructure=mapper.target.Select("cname=A and rnum=42,43,85")``
A :class:`RuntimeError` is raised if *substructure*
does not refer to the same underlying
:class:`ost.mol.EntityHandle` as :attr:`~target`.
:type substructure: :class:`ost.mol.EntityView`
:param model: Structure in which one wants to find representations for
*substructure*
:type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
:param topn: Max number of representations that are returned
:type topn: :class:`int`
:param inclusion_radius: Inclusion radius for lDDT
:type inclusion_radius: :class:`float`
:param thresholds: Thresholds for lDDT
:type thresholds: :class:`list` of :class:`float`
:param bb_only: Only consider backbone atoms in lDDT computation
:type bb_only: :class:`bool`
:param only_interchain: Only score interchain contacts in lDDT. Useful
if you want to identify interface patches.
:type only_interchain: :class:`bool`
:param chem_mapping_result: Pro param. The result of
:func:`~GetChemMapping` where you provided
*model*. If set, *model* parameter is not
used.
:type chem_mapping_result: :class:`tuple`
:param global_mapping: Pro param. Specify a global mapping result. This
fully defines the desired representation in the
model but extracts it and enriches it with all
the nice attributes of :class:`ReprResult`.
The target attribute in *global_mapping* must be
of the same entity as self.target and the model
attribute of *global_mapping* must be of the same
entity as *model*.
:type global_mapping: :class:`MappingResult`
:returns: :class:`list` of :class:`ReprResult`
"""
if topn < 1:
raise RuntimeError("topn must be >= 1")
if global_mapping is not None:
# ensure that this mapping is derived from the same structures
if global_mapping.target.handle.GetHashCode() != \
self.target.handle.GetHashCode():
raise RuntimeError("global_mapping.target must be the same "
"entity as self.target")
if global_mapping.model.handle.GetHashCode() != \
model.handle.GetHashCode():
raise RuntimeError("global_mapping.model must be the same "
"entity as model param")
# check whether substructure really is a subset of self.target
for r in substructure.residues:
ch_name = r.GetChain().GetName()
rnum = r.GetNumber()
target_r = self.target.FindResidue(ch_name, rnum)
if not target_r.IsValid():
raise RuntimeError(f"substructure has residue "
f"{r.GetQualifiedName()} which is not in "
f"self.target")
if target_r.handle.GetHashCode() != r.handle.GetHashCode():
raise RuntimeError(f"substructure has residue "
f"{r.GetQualifiedName()} which has an "
f"equivalent in self.target but it does "
f"not refer to the same underlying "
f"EntityHandle")
for a in r.atoms:
target_a = target_r.FindAtom(a.GetName())
if not target_a.IsValid():
raise RuntimeError(f"substructure has atom "
f"{a.GetQualifiedName()} which is not "
f"in self.target")
if a.handle.GetHashCode() != target_a.handle.GetHashCode():
raise RuntimeError(f"substructure has atom "
f"{a.GetQualifiedName()} which has an "
f"equivalent in self.target but it does "
f"not refer to the same underlying "
f"EntityHandle")
# check whether it contains either CA or C3'
ca = r.FindAtom("CA")
c3 = r.FindAtom("C3'") # FindAtom with prime in string is tested
# and works
if not ca.IsValid() and not c3.IsValid():
raise RuntimeError("All residues in substructure must contain "
"a backbone atom named CA or C3\'")
# perform mapping and alignments on full structures
if chem_mapping_result is None:
chem_mapping, chem_group_alns, mdl = self.GetChemMapping(model)
else:
chem_mapping, chem_group_alns, mdl = chem_mapping_result
ref_mdl_alns = _GetRefMdlAlns(self.chem_groups,
self.chem_group_alignments,
chem_mapping,
chem_group_alns)
# Get residue indices relative to full target chain
substructure_res_indices = dict()
for ch in substructure.chains:
full_ch = self.target.FindChain(ch.GetName())
idx = [full_ch.GetResidueIndex(r.GetNumber()) for r in ch.residues]
substructure_res_indices[ch.GetName()] = idx
# strip down variables to make them specific to substructure
# keep only chem_groups which are present in substructure
substructure_chem_groups = list()
substructure_chem_mapping = list()
chnames = set([ch.GetName() for ch in substructure.chains])
for chem_group, mapping in zip(self.chem_groups, chem_mapping):
substructure_chem_group = [ch for ch in chem_group if ch in chnames]
if len(substructure_chem_group) > 0:
substructure_chem_groups.append(substructure_chem_group)
substructure_chem_mapping.append(mapping)
# early stopping if no mdl chain can be mapped to substructure
n_mapped_mdl_chains = sum([len(m) for m in substructure_chem_mapping])
if n_mapped_mdl_chains == 0:
return list()
# strip the reference sequence in alignments to only contain
# sequence from substructure
substructure_ref_mdl_alns = dict()
mdl_views = dict()
for ch in mdl.chains:
mdl_views[ch.GetName()] = _CSel(mdl, [ch.GetName()])
for chem_group, mapping in zip(substructure_chem_groups,
substructure_chem_mapping):
for ref_ch in chem_group:
for mdl_ch in mapping:
full_aln = ref_mdl_alns[(ref_ch, mdl_ch)]
ref_seq = full_aln.GetSequence(0)
# the ref sequence is tricky... we start with a gap only
# sequence and only add olcs as defined by the residue
# indices that we extracted before...
tmp = ['-'] * len(full_aln)
for idx in substructure_res_indices[ref_ch]:
idx_in_seq = ref_seq.GetPos(idx)
tmp[idx_in_seq] = ref_seq[idx_in_seq]
ref_seq = seq.CreateSequence(ref_ch, ''.join(tmp))
ref_seq.AttachView(_CSel(substructure, [ref_ch]))
mdl_seq = full_aln.GetSequence(1)
mdl_seq = seq.CreateSequence(mdl_seq.GetName(),
mdl_seq.GetString())
mdl_seq.AttachView(mdl_views[mdl_ch])
aln = seq.CreateAlignment()
aln.AddSequence(ref_seq)
aln.AddSequence(mdl_seq)
substructure_ref_mdl_alns[(ref_ch, mdl_ch)] = aln
lddt_scorer = lddt.lDDTScorer(substructure,
inclusion_radius = inclusion_radius,
bb_only = bb_only)
scored_mappings = list()
if global_mapping:
# construct mapping of substructure from global mapping
flat_mapping = global_mapping.GetFlatMapping()
mapping = list()
for chem_group, chem_mapping in zip(substructure_chem_groups,
substructure_chem_mapping):
chem_group_mapping = list()
for ch in chem_group:
if ch in flat_mapping:
mdl_ch = flat_mapping[ch]
if mdl_ch in chem_mapping:
chem_group_mapping.append(mdl_ch)
else:
chem_group_mapping.append(None)
else:
chem_group_mapping.append(None)
mapping.append(chem_group_mapping)
mappings = [mapping]
else:
mappings = list(_ChainMappings(substructure_chem_groups,
substructure_chem_mapping,
self.n_max_naive))
for mapping in mappings:
# chain_mapping and alns as input for lDDT computation
lddt_chain_mapping = dict()
lddt_alns = dict()
n_res_aln = 0
for ref_chem_group, mdl_chem_group in zip(substructure_chem_groups,
mapping):
for ref_ch, mdl_ch in zip(ref_chem_group, mdl_chem_group):
# some mdl chains can be None
if mdl_ch is not None:
lddt_chain_mapping[mdl_ch] = ref_ch
aln = substructure_ref_mdl_alns[(ref_ch, mdl_ch)]
lddt_alns[mdl_ch] = aln
tmp = [int(c[0] != '-' and c[1] != '-') for c in aln]
n_res_aln += sum(tmp)
# don't compute lDDT if no single residue in mdl and ref is aligned
if n_res_aln == 0:
continue
lDDT, _ = lddt_scorer.lDDT(mdl, thresholds=thresholds,
chain_mapping=lddt_chain_mapping,
residue_mapping = lddt_alns,
check_resnames = False,
no_intrachain = only_interchain)
if lDDT is None:
ost.LogVerbose("No valid contacts in the reference")
lDDT = 0.0 # that means, that we have not a single valid contact
# in lDDT. For the code below to work, we just set it
# to a terrible score => 0.0
if len(scored_mappings) == 0:
scored_mappings.append((lDDT, mapping))
elif len(scored_mappings) < topn:
scored_mappings.append((lDDT, mapping))
scored_mappings.sort(reverse=True, key=lambda x: x[0])
elif lDDT > scored_mappings[-1][0]:
scored_mappings.append((lDDT, mapping))
scored_mappings.sort(reverse=True, key=lambda x: x[0])
scored_mappings = scored_mappings[:topn]
# finalize and return
results = list()
for scored_mapping in scored_mappings:
ref_view = substructure.handle.CreateEmptyView()
mdl_view = mdl.handle.CreateEmptyView()
for ref_ch_group, mdl_ch_group in zip(substructure_chem_groups,
scored_mapping[1]):
for ref_ch, mdl_ch in zip(ref_ch_group, mdl_ch_group):
if ref_ch is not None and mdl_ch is not None:
aln = substructure_ref_mdl_alns[(ref_ch, mdl_ch)]
for col in aln:
if col[0] != '-' and col[1] != '-':
ref_view.AddResidue(col.GetResidue(0),
mol.ViewAddFlag.INCLUDE_ALL)
mdl_view.AddResidue(col.GetResidue(1),
mol.ViewAddFlag.INCLUDE_ALL)
results.append(ReprResult(scored_mapping[0], substructure,
ref_view, mdl_view))
return results
def GetNMappings(self, model):
""" Returns number of possible mappings
:param model: Model with chains that are mapped onto
:attr:`chem_groups`
:type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
"""
chem_mapping, chem_group_alns, mdl = self.GetChemMapping(model)
return _NMappings(self.chem_groups, chem_mapping)
def ProcessStructure(self, ent):
""" Entity processing for chain mapping
* Selects view containing peptide and nucleotide residues which have
required backbone atoms present - for peptide residues thats
N, CA, C and CB (no CB for GLY), for nucleotide residues thats
O5', C5', C4', C3' and O3'.
* filters view by chain lengths, see *min_pep_length* and
*min_nuc_length* in constructor
* Extracts atom sequences for each chain in that view
* Attaches corresponding :class:`ost.mol.EntityView` to each sequence
* If residue number alignments are used, strictly increasing residue
numbers without insertion codes are ensured in each chain
:param ent: Entity to process
:type ent: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
:returns: Tuple with 3 elements: 1) :class:`ost.mol.EntityView`
containing peptide and nucleotide residues 2)
:class:`ost.seq.SequenceList` containing ATOMSEQ sequences
for each polypeptide chain in returned view, sequences have
:class:`ost.mol.EntityView` of according chains attached
3) same for polynucleotide chains
"""
view = ent.CreateEmptyView()
exp_pep_atoms = ["N", "CA", "C", "CB"]
exp_nuc_atoms = ["\"O5'\"", "\"C5'\"", "\"C4'\"", "\"C3'\"", "\"O3'\""]
pep_query = "peptide=true and aname=" + ','.join(exp_pep_atoms)
nuc_query = "nucleotide=true and aname=" + ','.join(exp_nuc_atoms)
pep_sel = ent.Select(pep_query)
for r in pep_sel.residues:
if len(r.atoms) == 4:
view.AddResidue(r.handle, mol.INCLUDE_ALL)
elif r.name == "GLY" and len(r.atoms) == 3:
atom_names = [a.GetName() for a in r.atoms]
if sorted(atom_names) == ["C", "CA", "N"]:
view.AddResidue(r.handle, mol.INCLUDE_ALL)
nuc_sel = ent.Select(nuc_query)
for r in nuc_sel.residues:
if len(r.atoms) == 5:
view.AddResidue(r.handle, mol.INCLUDE_ALL)
polypep_seqs = seq.CreateSequenceList()
polynuc_seqs = seq.CreateSequenceList()
if len(view.residues) == 0:
# no residues survived => return
return (view, polypep_seqs, polynuc_seqs)
for ch in view.chains:
n_res = len(ch.residues)
n_pep = sum([r.IsPeptideLinking() for r in ch.residues])
n_nuc = sum([r.IsNucleotideLinking() for r in ch.residues])
# guarantee that we have either pep or nuc (no mix of the two)
if n_pep > 0 and n_nuc > 0:
raise RuntimeError(f"Must not mix peptide and nucleotide linking "
f"residues in same chain ({ch.GetName()})")
if (n_pep + n_nuc) != n_res:
raise RuntimeError("All residues must either be peptide_linking "
"or nucleotide_linking")
# filter out short chains
if n_pep > 0 and n_pep < self.min_pep_length:
continue
if n_nuc > 0 and n_nuc < self.min_nuc_length:
continue
# the superfast residue number based alignment adds some
# restrictions on the numbers themselves:
# 1) no insertion codes 2) strictly increasing
if self.resnum_alignments:
# check if no insertion codes are present in residue numbers
ins_codes = [r.GetNumber().GetInsCode() for r in ch.residues]
if len(set(ins_codes)) != 1 or ins_codes[0] != '\0':
raise RuntimeError("Residue numbers in input structures must not "
"contain insertion codes")
# check if residue numbers are strictly increasing
nums = [r.GetNumber().GetNum() for r in ch.residues]
if not all(i < j for i, j in zip(nums, nums[1:])):
raise RuntimeError("Residue numbers in input structures must be "
"strictly increasing for each chain")
s = ''.join([r.one_letter_code for r in ch.residues])
s = seq.CreateSequence(ch.GetName(), s)
s.AttachView(_CSel(view, [ch.GetName()]))
if n_pep == n_res:
polypep_seqs.AddSequence(s)
elif n_nuc == n_res:
polynuc_seqs.AddSequence(s)
else:
raise RuntimeError("This shouldnt happen")
if len(polypep_seqs) == 0 and len(polynuc_seqs) == 0:
raise RuntimeError(f"No chain fulfilled minimum length requirement "
f"to be considered in chain mapping "
f"({self.min_pep_length} for peptide chains, "
f"{self.min_nuc_length} for nucleotide chains) "
f"- mapping failed")
# select for chains for which we actually extracted the sequence
chain_names = [s.GetAttachedView().chains[0].name for s in polypep_seqs]
chain_names += [s.GetAttachedView().chains[0].name for s in polynuc_seqs]
view = _CSel(view, chain_names)
return (view, polypep_seqs, polynuc_seqs)
def Align(self, s1, s2, stype):
""" Access to internal sequence alignment functionality
Alignment parameterization is setup at ChainMapper construction
:param s1: First sequence to align - must have view attached in case
of resnum_alignments
:type s1: :class:`ost.seq.SequenceHandle`
:param s2: Second sequence to align - must have view attached in case
of resnum_alignments
:type s2: :class:`ost.seq.SequenceHandle`
:param stype: Type of sequences to align, must be in
[:class:`ost.mol.ChemType.AMINOACIDS`,
:class:`ost.mol.ChemType.NUCLEOTIDES`]
:returns: Pairwise alignment of s1 and s2
"""
if stype not in [mol.ChemType.AMINOACIDS, mol.ChemType.NUCLEOTIDES]:
raise RuntimeError("stype must be ost.mol.ChemType.AMINOACIDS or "
"ost.mol.ChemType.NUCLEOTIDES")
return self.aligner.Align(s1, s2, chem_type = stype)
# INTERNAL HELPERS
##################
class _Aligner:
def __init__(self, pep_subst_mat = seq.alg.BLOSUM62, pep_gap_open = -5,
pep_gap_ext = -2, nuc_subst_mat = seq.alg.NUC44,
nuc_gap_open = -4, nuc_gap_ext = -4, resnum_aln = False):
""" Helper class to compute alignments
Sets default values for substitution matrix, gap open and gap extension
penalties. They are only used in default mode (Needleman-Wunsch aln).
If *resnum_aln* is True, only residue numbers of views that are attached
to input sequences are considered.
"""
self.pep_subst_mat = pep_subst_mat
self.pep_gap_open = pep_gap_open
self.pep_gap_ext = pep_gap_ext
self.nuc_subst_mat = nuc_subst_mat
self.nuc_gap_open = nuc_gap_open
self.nuc_gap_ext = nuc_gap_ext
self.resnum_aln = resnum_aln
def Align(self, s1, s2, chem_type=None):
if self.resnum_aln:
return self.ResNumAlign(s1, s2)
else:
if chem_type is None:
raise RuntimeError("Must specify chem_type for NW alignment")
return self.NWAlign(s1, s2, chem_type)
def NWAlign(self, s1, s2, chem_type):
""" Returns pairwise alignment using Needleman-Wunsch algorithm
:param s1: First sequence to align
:type s1: :class:`ost.seq.SequenceHandle`
:param s2: Second sequence to align
:type s2: :class:`ost.seq.SequenceHandle`
:param chem_type: Must be in [:class:`ost.mol.ChemType.AMINOACIDS`,
:class:`ost.mol.ChemType.NUCLEOTIDES`], determines
substitution matrix and gap open/extension penalties
:type chem_type: :class:`ost.mol.ChemType`
:returns: Alignment with s1 as first and s2 as second sequence
"""
if chem_type == mol.ChemType.AMINOACIDS:
return seq.alg.GlobalAlign(s1, s2, self.pep_subst_mat,
gap_open=self.pep_gap_open,
gap_ext=self.pep_gap_ext)[0]
elif chem_type == mol.ChemType.NUCLEOTIDES:
return seq.alg.GlobalAlign(s1, s2, self.nuc_subst_mat,
gap_open=self.nuc_gap_open,
gap_ext=self.nuc_gap_ext)[0]
else:
raise RuntimeError("Invalid ChemType")
return aln
def ResNumAlign(self, s1, s2):
""" Returns pairwise alignment using residue numbers of attached views
Assumes that there are no insertion codes (alignment only on numerical
component) and that resnums are strictly increasing (fast min/max
identification). These requirements are assured if a structure has been
processed by :class:`ChainMapper`.
:param s1: First sequence to align, must have :class:`ost.mol.EntityView`
attached
:type s1: :class:`ost.seq.SequenceHandle`
:param s2: Second sequence to align, must have :class:`ost.mol.EntityView`
attached
:type s2: :class:`ost.seq.SequenceHandle`
"""
assert(s1.HasAttachedView())
assert(s2.HasAttachedView())
v1 = s1.GetAttachedView()
rnums1 = [r.GetNumber().GetNum() for r in v1.residues]
v2 = s2.GetAttachedView()
rnums2 = [r.GetNumber().GetNum() for r in v2.residues]
min_num = min(rnums1[0], rnums2[0])
max_num = max(rnums1[-1], rnums2[-1])
aln_length = max_num - min_num + 1
aln_s1 = ['-'] * aln_length
for r, rnum in zip(v1.residues, rnums1):
aln_s1[rnum-min_num] = r.one_letter_code
aln_s2 = ['-'] * aln_length
for r, rnum in zip(v2.residues, rnums2):
aln_s2[rnum-min_num] = r.one_letter_code
aln = seq.CreateAlignment()
aln.AddSequence(seq.CreateSequence(s1.GetName(), ''.join(aln_s1)))
aln.AddSequence(seq.CreateSequence(s2.GetName(), ''.join(aln_s2)))
return aln
def _GetAlnPropsTwo(aln):
"""Returns basic properties of *aln* version two...
:param aln: Alignment to compute properties
:type aln: :class:`seq.AlignmentHandle`
:returns: Tuple with 2 elements. 1) sequence identify in range [0, 100]
considering aligned columns 2) Fraction of non-gap characters
in first sequence that are covered by non-gap characters in
second sequence.
"""
assert(aln.GetCount() == 2)
n_tot = sum([1 for col in aln if col[0] != '-'])
n_aligned = sum([1 for col in aln if (col[0] != '-' and col[1] != '-')])
return (seq.alg.SequenceIdentity(aln), float(n_aligned)/n_tot)
def _GetAlnPropsOne(aln):
"""Returns basic properties of *aln* version one...
:param aln: Alignment to compute properties
:type aln: :class:`seq.AlignmentHandle`
:returns: Tuple with 3 elements. 1) sequence identify in range [0, 100]
considering aligned columns 2) Fraction of gaps between
first and last aligned column in s1 3) same for s2.
"""
assert(aln.GetCount() == 2)
n_gaps_1 = str(aln.GetSequence(0)).strip('-').count('-')
n_gaps_2 = str(aln.GetSequence(1)).strip('-').count('-')
gap_frac_1 = float(n_gaps_1)/len(aln.GetSequence(0).GetGaplessString())
gap_frac_2 = float(n_gaps_2)/len(aln.GetSequence(1).GetGaplessString())
return (seq.alg.SequenceIdentity(aln), gap_frac_1, gap_frac_2)
def _GetChemGroupAlignments(pep_seqs, nuc_seqs, aligner, pep_seqid_thr=95.,
pep_gap_thr=0.1, nuc_seqid_thr=95.,
nuc_gap_thr=0.1):
"""Returns alignments with groups of chemically equivalent chains
:param pep_seqs: List of polypeptide sequences
:type pep_seqs: :class:`seq.SequenceList`
:param nuc_seqs: List of polynucleotide sequences
:type nuc_seqs: :class:`seq.SequenceList`
:param aligner: Helper class to generate pairwise alignments
:type aligner: :class:`_Aligner`
:param pep_seqid_thr: Threshold used to decide when two peptide chains are
identical. 95 percent tolerates the few mutations
crystallographers like to do.
:type pep_seqid_thr: :class:`float`
:param pep_gap_thr: Additional threshold to avoid gappy alignments with high
seqid. The reason for not just normalizing seqid by the
longer sequence is that one sequence might be a perfect
subsequence of the other but only cover half of it. This
threshold checks for a maximum allowed fraction of gaps
in any of the two sequences after stripping terminal gaps.
:type pep_gap_thr: :class:`float`
:param nuc_seqid_thr: Nucleotide equivalent of *pep_seqid_thr*
:type nuc_seqid_thr: :class:`float`
:param nuc_gap_thr: Nucleotide equivalent of *nuc_gap_thr*
:type nuc_gap_thr: :class:`float`
:returns: Tuple with first element being an AlignmentList. Each alignment
represents a group of chemically equivalent chains and the first
sequence is the longest. Second element is a list of equivalent
length specifying the types of the groups. List elements are in
[:class:`ost.ChemType.AMINOACIDS`,
:class:`ost.ChemType.NUCLEOTIDES`]
"""
pep_groups = _GroupSequences(pep_seqs, pep_seqid_thr, pep_gap_thr, aligner,
mol.ChemType.AMINOACIDS)
nuc_groups = _GroupSequences(nuc_seqs, nuc_seqid_thr, nuc_gap_thr, aligner,
mol.ChemType.NUCLEOTIDES)
group_types = [mol.ChemType.AMINOACIDS] * len(pep_groups)
group_types += [mol.ChemType.NUCLEOTIDES] * len(nuc_groups)
groups = pep_groups
groups.extend(nuc_groups)
return (groups, group_types)
def _GroupSequences(seqs, seqid_thr, gap_thr, aligner, chem_type):
"""Get list of alignments representing groups of equivalent sequences
:param seqid_thr: Threshold used to decide when two chains are identical.
:type seqid_thr: :class:`float`
:param gap_thr: Additional threshold to avoid gappy alignments with high
seqid. The reason for not just normalizing seqid by the
longer sequence is that one sequence might be a perfect
subsequence of the other but only cover half of it. This
threshold checks for a maximum allowed fraction of gaps
in any of the two sequences after stripping terminal gaps.
:type gap_thr: :class:`float`
:param aligner: Helper class to generate pairwise alignments
:type aligner: :class:`_Aligner`
:param chem_type: ChemType of seqs which is passed to *aligner*, must be in
[:class:`ost.mol.ChemType.AMINOACIDS`,
:class:`ost.mol.ChemType.NUCLEOTIDES`]
:type chem_type: :class:`ost.mol.ChemType`
:returns: A list of alignments, one alignment for each group
with longest sequence (reference) as first sequence.
:rtype: :class:`ost.seq.AlignmentList`
"""
groups = list()
for s_idx in range(len(seqs)):
matching_group = None
for g_idx in range(len(groups)):
for g_s_idx in range(len(groups[g_idx])):
aln = aligner.Align(seqs[s_idx], seqs[groups[g_idx][g_s_idx]],
chem_type)
sid, frac_i, frac_j = _GetAlnPropsOne(aln)
if sid >= seqid_thr and frac_i < gap_thr and frac_j < gap_thr:
matching_group = g_idx
break
if matching_group is not None:
break
if matching_group is None:
groups.append([s_idx])
else:
groups[matching_group].append(s_idx)
# sort based on sequence length
sorted_groups = list()
for g in groups:
if len(g) > 1:
tmp = sorted([[len(seqs[i]), i] for i in g], reverse=True)
sorted_groups.append([x[1] for x in tmp])
else:
sorted_groups.append(g)
# translate from indices back to sequences and directly generate alignments
# of the groups with the longest (first) sequence as reference
aln_list = seq.AlignmentList()
for g in sorted_groups:
if len(g) == 1:
# aln with one single sequence
aln_list.append(seq.CreateAlignment(seqs[g[0]]))
else:
# obtain pairwise aln of first sequence (reference) to all others
alns = seq.AlignmentList()
i = g[0]
for j in g[1:]:
alns.append(aligner.Align(seqs[i], seqs[j], chem_type))
# and merge
aln_list.append(seq.alg.MergePairwiseAlignments(alns, seqs[i]))
# transfer attached views
seq_dict = {s.GetName(): s for s in seqs}
for aln_idx in range(len(aln_list)):
for aln_s_idx in range(aln_list[aln_idx].GetCount()):
s_name = aln_list[aln_idx].GetSequence(aln_s_idx).GetName()
s = seq_dict[s_name]
aln_list[aln_idx].AttachView(aln_s_idx, s.GetAttachedView())
return aln_list
def _MapSequence(ref_seqs, ref_types, s, s_type, aligner):
"""Tries top map *s* onto any of the sequences in *ref_seqs*
Computes alignments of *s* to each of the reference sequences of equal type
and sorts them by seqid*fraction_covered (seqid: sequence identity of
aligned columns in alignment, fraction_covered: Fraction of non-gap
characters in reference sequence that are covered by non-gap characters in
*s*). Best scoring mapping is returned.
:param ref_seqs: Reference sequences
:type ref_seqs: :class:`ost.seq.SequenceList`
:param ref_types: Types of reference sequences, e.g.
ost.mol.ChemType.AminoAcids
:type ref_types: :class:`list` of :class:`ost.mol.ChemType`
:param s: Sequence to map
:type s: :class:`ost.seq.SequenceHandle`
:param s_type: Type of *s*, only try mapping to sequences in *ref_seqs*
with equal type as defined in *ref_types*
:param aligner: Helper class to generate pairwise alignments
:type aligner: :class:`_Aligner`
:returns: Tuple with two elements. 1) index of sequence in *ref_seqs* to
which *s* can be mapped 2) Pairwise sequence alignment with
sequence from *ref_seqs* as first sequence. Both elements are
None if no mapping can be found.
:raises: :class:`RuntimeError` if mapping is ambiguous, i.e. *s*
successfully maps to more than one sequence in *ref_seqs*
"""
scored_alns = list()
for ref_idx, ref_seq in enumerate(ref_seqs):
if ref_types[ref_idx] == s_type:
aln = aligner.Align(ref_seq, s, s_type)
seqid, fraction_covered = _GetAlnPropsTwo(aln)
score = seqid * fraction_covered
scored_alns.append((score, ref_idx, aln))
if len(scored_alns) == 0:
return (None, None) # no mapping possible...
scored_alns = sorted(scored_alns, key=lambda x: x[0], reverse=True)
return (scored_alns[0][1], scored_alns[0][2])
def _GetRefMdlAlns(ref_chem_groups, ref_chem_group_msas, mdl_chem_groups,
mdl_chem_group_alns, pairs=None):
""" Get all possible ref/mdl chain alignments given chem group mapping
:param ref_chem_groups: :attr:`ChainMapper.chem_groups`
:type ref_chem_groups: :class:`list` of :class:`list` of :class:`str`
:param ref_chem_group_msas: :attr:`ChainMapper.chem_group_alignments`
:type ref_chem_group_msas: :class:`ost.seq.AlignmentList`
:param mdl_chem_groups: Groups of model chains that are mapped to
*ref_chem_groups*. Return value of
:func:`ChainMapper.GetChemMapping`.
:type mdl_chem_groups: :class:`list` of :class:`list` of :class:`str`
:param mdl_chem_group_alns: A pairwise sequence alignment for every chain
in *mdl_chem_groups* that aligns these sequences
to the respective reference sequence.
Return values of
:func:`ChainMapper.GetChemMapping`.
:type mdl_chem_group_alns: :class:`list` of :class:`ost.seq.AlignmentList`
:param pairs: Pro param - restrict return dict to specified pairs. A set of
tuples in form (<trg_ch>, <mdl_ch>)
:type pairs: :class:`set`
:returns: A dictionary holding all possible ref/mdl chain alignments. Keys
in that dictionary are tuples of the form (ref_ch, mdl_ch) and
values are the respective pairwise alignments with first sequence
being from ref, the second from mdl.
"""
# alignment of each model chain to chem_group reference sequence
mdl_alns = dict()
for alns in mdl_chem_group_alns:
for aln in alns:
mdl_chain_name = aln.GetSequence(1).GetName()
mdl_alns[mdl_chain_name] = aln
# generate all alignments between ref/mdl chain atomseqs that we will
# ever observe
ref_mdl_alns = dict()
for ref_chains, mdl_chains, ref_aln in zip(ref_chem_groups, mdl_chem_groups,
ref_chem_group_msas):
for ref_ch in ref_chains:
for mdl_ch in mdl_chains:
if pairs is not None and (ref_ch, mdl_ch) not in pairs:
continue
# obtain alignments of mdl and ref chains towards chem
# group ref sequence and merge them
aln_list = seq.AlignmentList()
# do ref aln
s1 = ref_aln.GetSequence(0)
s2 = ref_aln.GetSequence(ref_chains.index(ref_ch))
aln_list.append(seq.CreateAlignment(s1, s2))
# do mdl aln
aln_list.append(mdl_alns[mdl_ch])
# merge
ref_seq = seq.CreateSequence(s1.GetName(),
s1.GetGaplessString())
merged_aln = seq.alg.MergePairwiseAlignments(aln_list,
ref_seq)
# merged_aln:
# seq1: ref seq of chem group
# seq2: seq of ref chain
# seq3: seq of mdl chain
# => we need the alignment between seq2 and seq3
s2 = merged_aln.GetSequence(1)
s3 = merged_aln.GetSequence(2)
# cut leading and trailing gap columns
a = 0 # number of leading gap columns
for idx in range(len(s2)):
if s2[idx] != '-' or s3[idx] != '-':
break
a += 1
b = 0 # number of trailing gap columns
for idx in reversed(range(len(s2))):
if s2[idx] != '-' or s3[idx] != '-':
break
b += 1
s2 = seq.CreateSequence(s2.GetName(), s2[a: len(s2)-b])
s3 = seq.CreateSequence(s3.GetName(), s3[a: len(s3)-b])
ref_mdl_alns[(ref_ch, mdl_ch)] = seq.CreateAlignment(s2, s3)
return ref_mdl_alns
def _CheckOneToOneMapping(ref_chains, mdl_chains):
""" Checks whether we already have a perfect one to one mapping
That means each list in *ref_chains* has exactly one element and each
list in *mdl_chains* has either one element (it's mapped) or is empty
(ref chain has no mapped mdl chain). Returns None if no such mapping
can be found.
:param ref_chains: corresponds to :attr:`ChainMapper.chem_groups`
:type ref_chains: :class:`list` of :class:`list` of :class:`str`
:param mdl_chains: mdl chains mapped to chem groups in *ref_chains*, i.e.
the return value of :func:`ChainMapper.GetChemMapping`
:type mdl_chains: class:`list` of :class:`list` of :class:`str`
:returns: A :class:`list` of :class:`list` if a one to one mapping is found,
None otherwise
"""
only_one_to_one = True
one_to_one = list()
for ref, mdl in zip(ref_chains, mdl_chains):
if len(ref) == 1 and len(mdl) == 1:
one_to_one.append(mdl)
elif len(ref) == 1 and len(mdl) == 0:
one_to_one.append([None])
else:
only_one_to_one = False
break
if only_one_to_one:
return one_to_one
else:
return None
class _lDDTDecomposer:
def __init__(self, ref, mdl, ref_mdl_alns, inclusion_radius = 15.0,
thresholds = [0.5, 1.0, 2.0, 4.0]):
""" Compute backbone only lDDT scores for ref/mdl
Uses the pairwise decomposable property of backbone only lDDT and
implements a caching mechanism to efficiently enumerate different
chain mappings.
"""
self.ref = ref
self.mdl = mdl
self.ref_mdl_alns = ref_mdl_alns
self.inclusion_radius = inclusion_radius
self.thresholds = thresholds
# keep track of single chains and interfaces in ref
self.ref_chains = list() # e.g. ['A', 'B', 'C']
self.ref_interfaces = list() # e.g. [('A', 'B'), ('A', 'C')]
# holds lDDT scorer for each chain in ref
# key: chain name, value: scorer
self.single_chain_scorer = dict()
# cache for single chain conserved contacts
# key: tuple (ref_ch, mdl_ch) value: number of conserved contacts
self.single_chain_cache = dict()
# holds lDDT scorer for each pairwise interface in target
# key: tuple (ref_ch1, ref_ch2), value: scorer
self.interface_scorer = dict()
# cache for interface conserved contacts
# key: tuple of tuple ((ref_ch1, ref_ch2),((mdl_ch1, mdl_ch2))
# value: number of conserved contacts
self.interface_cache = dict()
self.n = 0
self._SetupScorer()
def _SetupScorer(self):
for ch in self.ref.chains:
# Select everything close to that chain
query = f"{self.inclusion_radius} <> "
query += f"[cname={mol.QueryQuoteName(ch.GetName())}] "
query += f"and cname!={mol.QueryQuoteName(ch.GetName())}"
for close_ch in self.ref.Select(query).chains:
k1 = (ch.GetName(), close_ch.GetName())
k2 = (close_ch.GetName(), ch.GetName())
if k1 not in self.interface_scorer and \
k2 not in self.interface_scorer:
dimer_ref = _CSel(self.ref, [k1[0], k1[1]])
s = lddt.lDDTScorer(dimer_ref, bb_only=True)
self.interface_scorer[k1] = s
self.interface_scorer[k2] = s
self.n += self.interface_scorer[k1].n_distances_ic
self.ref_interfaces.append(k1)
# single chain scorer are actually interface scorers to save
# some distance calculations
if ch.GetName() not in self.single_chain_scorer:
self.single_chain_scorer[ch.GetName()] = s
self.n += s.GetNChainContacts(ch.GetName(),
no_interchain=True)
self.ref_chains.append(ch.GetName())
if close_ch.GetName() not in self.single_chain_scorer:
self.single_chain_scorer[close_ch.GetName()] = s
self.n += s.GetNChainContacts(close_ch.GetName(),
no_interchain=True)
self.ref_chains.append(close_ch.GetName())
# add any missing single chain scorer
for ch in self.ref.chains:
if ch.GetName() not in self.single_chain_scorer:
single_chain_ref = _CSel(self.ref, [ch.GetName()])
self.single_chain_scorer[ch.GetName()] = \
lddt.lDDTScorer(single_chain_ref, bb_only = True)
self.n += self.single_chain_scorer[ch.GetName()].n_distances
self.ref_chains.append(ch.GetName())
def lDDT(self, ref_chain_groups, mdl_chain_groups):
flat_map = dict()
for ref_chains, mdl_chains in zip(ref_chain_groups, mdl_chain_groups):
for ref_ch, mdl_ch in zip(ref_chains, mdl_chains):
flat_map[ref_ch] = mdl_ch
return self.lDDTFromFlatMap(flat_map)
def lDDTFromFlatMap(self, flat_map):
conserved = 0
# do single chain scores
for ref_ch in self.ref_chains:
if ref_ch in flat_map and flat_map[ref_ch] is not None:
conserved += self.SCCounts(ref_ch, flat_map[ref_ch])
# do interfaces
for ref_ch1, ref_ch2 in self.ref_interfaces:
if ref_ch1 in flat_map and ref_ch2 in flat_map:
mdl_ch1 = flat_map[ref_ch1]
mdl_ch2 = flat_map[ref_ch2]
if mdl_ch1 is not None and mdl_ch2 is not None:
conserved += self.IntCounts(ref_ch1, ref_ch2, mdl_ch1,
mdl_ch2)
return conserved / (len(self.thresholds) * self.n)
def SCCounts(self, ref_ch, mdl_ch):
if not (ref_ch, mdl_ch) in self.single_chain_cache:
alns = dict()
alns[mdl_ch] = self.ref_mdl_alns[(ref_ch, mdl_ch)]
mdl_sel = _CSel(self.mdl, [mdl_ch])
s = self.single_chain_scorer[ref_ch]
_,_,_,conserved,_,_,_ = s.lDDT(mdl_sel,
residue_mapping=alns,
return_dist_test=True,
no_interchain=True,
chain_mapping={mdl_ch: ref_ch},
check_resnames=False)
self.single_chain_cache[(ref_ch, mdl_ch)] = conserved
return self.single_chain_cache[(ref_ch, mdl_ch)]
def IntCounts(self, ref_ch1, ref_ch2, mdl_ch1, mdl_ch2):
k1 = ((ref_ch1, ref_ch2),(mdl_ch1, mdl_ch2))
k2 = ((ref_ch2, ref_ch1),(mdl_ch2, mdl_ch1))
if k1 not in self.interface_cache and k2 not in self.interface_cache:
alns = dict()
alns[mdl_ch1] = self.ref_mdl_alns[(ref_ch1, mdl_ch1)]
alns[mdl_ch2] = self.ref_mdl_alns[(ref_ch2, mdl_ch2)]
mdl_sel = _CSel(self.mdl, [mdl_ch1, mdl_ch2])
s = self.interface_scorer[(ref_ch1, ref_ch2)]
_,_,_,conserved,_,_,_ = s.lDDT(mdl_sel,
residue_mapping=alns,
return_dist_test=True,
no_intrachain=True,
chain_mapping={mdl_ch1: ref_ch1,
mdl_ch2: ref_ch2},
check_resnames=False)
self.interface_cache[k1] = conserved
self.interface_cache[k2] = conserved
return self.interface_cache[k1]
class _lDDTGreedySearcher(_lDDTDecomposer):
def __init__(self, ref, mdl, ref_chem_groups, mdl_chem_groups,
ref_mdl_alns, inclusion_radius = 15.0,
thresholds = [0.5, 1.0, 2.0, 4.0],
steep_opt_rate = None):
""" Greedy extension of already existing but incomplete chain mappings
"""
super().__init__(ref, mdl, ref_mdl_alns,
inclusion_radius = inclusion_radius,
thresholds = thresholds)
self.steep_opt_rate = steep_opt_rate
self.neighbors = {k: set() for k in self.ref_chains}
for k in self.interface_scorer.keys():
self.neighbors[k[0]].add(k[1])
self.neighbors[k[1]].add(k[0])
assert(len(ref_chem_groups) == len(mdl_chem_groups))
self.ref_chem_groups = ref_chem_groups
self.mdl_chem_groups = mdl_chem_groups
self.ref_ch_group_mapper = dict()
self.mdl_ch_group_mapper = dict()
for g_idx, (ref_g, mdl_g) in enumerate(zip(ref_chem_groups,
mdl_chem_groups)):
for ch in ref_g:
self.ref_ch_group_mapper[ch] = g_idx
for ch in mdl_g:
self.mdl_ch_group_mapper[ch] = g_idx
# keep track of mdl chains that potentially give lDDT contributions,
# i.e. they have locations within inclusion_radius + max(thresholds)
self.mdl_neighbors = dict()
d = self.inclusion_radius + max(self.thresholds)
for ch in self.mdl.chains:
ch_name = ch.GetName()
self.mdl_neighbors[ch_name] = set()
query = f"{d} <> [cname={mol.QueryQuoteName(ch_name)}]"
query += f" and cname !={mol.QueryQuoteName(ch_name)}"
for close_ch in self.mdl.Select(query).chains:
self.mdl_neighbors[ch_name].add(close_ch.GetName())
def ExtendMapping(self, mapping, max_ext = None):
if len(mapping) == 0:
raise RuntimError("Mapping must contain a starting point")
# Ref chains onto which we can map. The algorithm starts with a mapping
# on ref_ch. From there we can start to expand to connected neighbors.
# All neighbors that we can reach from the already mapped chains are
# stored in this set which will be updated during runtime
map_targets = set()
for ref_ch in mapping.keys():
map_targets.update(self.neighbors[ref_ch])
# remove the already mapped chains
for ref_ch in mapping.keys():
map_targets.discard(ref_ch)
if len(map_targets) == 0:
return mapping # nothing to extend
# keep track of what model chains are not yet mapped for each chem group
free_mdl_chains = list()
for chem_group in self.mdl_chem_groups:
tmp = [x for x in chem_group if x not in mapping.values()]
free_mdl_chains.append(set(tmp))
# keep track of what ref chains got a mapping
newly_mapped_ref_chains = list()
something_happened = True
while something_happened:
something_happened=False
if self.steep_opt_rate is not None:
n_chains = len(newly_mapped_ref_chains)
if n_chains > 0 and n_chains % self.steep_opt_rate == 0:
mapping = self._SteepOpt(mapping, newly_mapped_ref_chains)
if max_ext is not None and len(newly_mapped_ref_chains) >= max_ext:
break
max_n = 0
max_mapping = None
for ref_ch in map_targets:
chem_group_idx = self.ref_ch_group_mapper[ref_ch]
for mdl_ch in free_mdl_chains[chem_group_idx]:
# single chain score
n_single = self.SCCounts(ref_ch, mdl_ch)
# scores towards neighbors that are already mapped
n_inter = 0
for neighbor in self.neighbors[ref_ch]:
if neighbor in mapping and mapping[neighbor] in \
self.mdl_neighbors[mdl_ch]:
n_inter += self.IntCounts(ref_ch, neighbor, mdl_ch,
mapping[neighbor])
n = n_single + n_inter
if n_inter > 0 and n > max_n:
# Only accept a new solution if its actually connected
# i.e. n_inter > 0. Otherwise we could just map a big
# fat mdl chain sitting somewhere in Nirvana
max_n = n
max_mapping = (ref_ch, mdl_ch)
if max_n > 0:
something_happened = True
# assign new found mapping
mapping[max_mapping[0]] = max_mapping[1]
# add all neighboring chains to map targets as they are now
# reachable
for neighbor in self.neighbors[max_mapping[0]]:
if neighbor not in mapping:
map_targets.add(neighbor)
# remove the ref chain from map targets
map_targets.remove(max_mapping[0])
# remove the mdl chain from free_mdl_chains - its taken...
chem_group_idx = self.ref_ch_group_mapper[max_mapping[0]]
free_mdl_chains[chem_group_idx].remove(max_mapping[1])
# keep track of what ref chains got a mapping
newly_mapped_ref_chains.append(max_mapping[0])
return mapping
def _SteepOpt(self, mapping, chains_to_optimize=None):
# just optimize ALL ref chains if nothing specified
if chains_to_optimize is None:
chains_to_optimize = mapping.keys()
# make sure that we only have ref chains which are actually mapped
ref_chains = [x for x in chains_to_optimize if mapping[x] is not None]
# group ref chains to be optimized into chem groups
tmp = dict()
for ch in ref_chains:
chem_group_idx = self.ref_ch_group_mapper[ch]
if chem_group_idx in tmp:
tmp[chem_group_idx].append(ch)
else:
tmp[chem_group_idx] = [ch]
chem_groups = list(tmp.values())
# try all possible mapping swaps. Swaps that improve the score are
# immediately accepted and we start all over again
current_lddt = self.lDDTFromFlatMap(mapping)
something_happened = True
while something_happened:
something_happened = False
for chem_group in chem_groups:
if something_happened:
break
for ch1, ch2 in itertools.combinations(chem_group, 2):
swapped_mapping = dict(mapping)
swapped_mapping[ch1] = mapping[ch2]
swapped_mapping[ch2] = mapping[ch1]
score = self.lDDTFromFlatMap(swapped_mapping)
if score > current_lddt:
something_happened = True
mapping = swapped_mapping
current_lddt = score
break
return mapping
def _lDDTNaive(trg, mdl, inclusion_radius, thresholds, chem_groups,
chem_mapping, ref_mdl_alns, n_max_naive):
""" Naively iterates all possible chain mappings and returns the best
"""
best_mapping = None
best_lddt = -1.0
# Benchmarks on homo-oligomers indicate that full blown lDDT
# computation is faster up to tetramers => 4!=24 possible mappings.
# For stuff bigger than that, the decomposer approach should be used
if _NMappingsWithin(chem_groups, chem_mapping, 24):
# Setup scoring
lddt_scorer = lddt.lDDTScorer(trg, bb_only = True)
for mapping in _ChainMappings(chem_groups, chem_mapping, n_max_naive):
# chain_mapping and alns as input for lDDT computation
lddt_chain_mapping = dict()
lddt_alns = dict()
for ref_chem_group, mdl_chem_group in zip(chem_groups, mapping):
for ref_ch, mdl_ch in zip(ref_chem_group, mdl_chem_group):
# some mdl chains can be None
if mdl_ch is not None:
lddt_chain_mapping[mdl_ch] = ref_ch
lddt_alns[mdl_ch] = ref_mdl_alns[(ref_ch, mdl_ch)]
lDDT, _ = lddt_scorer.lDDT(mdl, thresholds=thresholds,
chain_mapping=lddt_chain_mapping,
residue_mapping = lddt_alns,
check_resnames = False)
if lDDT > best_lddt:
best_mapping = mapping
best_lddt = lDDT
else:
# Setup scoring
lddt_scorer = _lDDTDecomposer(trg, mdl, ref_mdl_alns,
inclusion_radius=inclusion_radius,
thresholds = thresholds)
for mapping in _ChainMappings(chem_groups, chem_mapping, n_max_naive):
lDDT = lddt_scorer.lDDT(chem_groups, mapping)
if lDDT > best_lddt:
best_mapping = mapping
best_lddt = lDDT
return (best_mapping, best_lddt)
def _GetSeeds(ref_chem_groups, mdl_chem_groups, mapped_ref_chains = set(),
mapped_mdl_chains = set()):
seeds = list()
for ref_chains, mdl_chains in zip(ref_chem_groups,
mdl_chem_groups):
for ref_ch in ref_chains:
if ref_ch not in mapped_ref_chains:
for mdl_ch in mdl_chains:
if mdl_ch not in mapped_mdl_chains:
seeds.append((ref_ch, mdl_ch))
return seeds
def _lDDTGreedyFast(the_greed):
something_happened = True
mapping = dict()
while something_happened:
something_happened = False
seeds = _GetSeeds(the_greed.ref_chem_groups,
the_greed.mdl_chem_groups,
mapped_ref_chains = set(mapping.keys()),
mapped_mdl_chains = set(mapping.values()))
# search for best scoring starting point
n_best = 0
best_seed = None
for seed in seeds:
n = the_greed.SCCounts(seed[0], seed[1])
if n > n_best:
n_best = n
best_seed = seed
if n_best == 0:
break # no proper seed found anymore...
# add seed to mapping and start the greed
mapping[best_seed[0]] = best_seed[1]
mapping = the_greed.ExtendMapping(mapping)
something_happened = True
# translate mapping format and return
final_mapping = list()
for ref_chains in the_greed.ref_chem_groups:
mapped_mdl_chains = list()
for ref_ch in ref_chains:
if ref_ch in mapping:
mapped_mdl_chains.append(mapping[ref_ch])
else:
mapped_mdl_chains.append(None)
final_mapping.append(mapped_mdl_chains)
return final_mapping
def _lDDTGreedyFull(the_greed):
""" Uses each reference chain as starting point for expansion
"""
seeds = _GetSeeds(the_greed.ref_chem_groups, the_greed.mdl_chem_groups)
best_overall_score = -1.0
best_overall_mapping = dict()
for seed in seeds:
# do initial extension
mapping = the_greed.ExtendMapping({seed[0]: seed[1]})
# repeat the process until we have a full mapping
something_happened = True
while something_happened:
something_happened = False
remnant_seeds = _GetSeeds(the_greed.ref_chem_groups,
the_greed.mdl_chem_groups,
mapped_ref_chains = set(mapping.keys()),
mapped_mdl_chains = set(mapping.values()))
if len(remnant_seeds) > 0:
# still more mapping to be done
best_score = -1.0
best_mapping = None
for remnant_seed in remnant_seeds:
tmp_mapping = dict(mapping)
tmp_mapping[remnant_seed[0]] = remnant_seed[1]
tmp_mapping = the_greed.ExtendMapping(tmp_mapping)
score = the_greed.lDDTFromFlatMap(tmp_mapping)
if score > best_score:
best_score = score
best_mapping = tmp_mapping
if best_mapping is not None:
something_happened = True
mapping = best_mapping
score = the_greed.lDDTFromFlatMap(mapping)
if score > best_overall_score:
best_overall_score = score
best_overall_mapping = mapping
mapping = best_overall_mapping
# translate mapping format and return
final_mapping = list()
for ref_chains in the_greed.ref_chem_groups:
mapped_mdl_chains = list()
for ref_ch in ref_chains:
if ref_ch in mapping:
mapped_mdl_chains.append(mapping[ref_ch])
else:
mapped_mdl_chains.append(None)
final_mapping.append(mapped_mdl_chains)
return final_mapping
def _lDDTGreedyBlock(the_greed, seed_size, blocks_per_chem_group):
""" try multiple seeds, i.e. try all ref/mdl chain combinations within the
respective chem groups and compute single chain lDDTs. The
*blocks_per_chem_group* best scoring ones are extend by *seed_size* chains
and the best scoring one is exhaustively extended.
"""
if seed_size is None or seed_size < 1:
raise RuntimeError(f"seed_size must be an int >= 1 (got {seed_size})")
if blocks_per_chem_group is None or blocks_per_chem_group < 1:
raise RuntimeError(f"blocks_per_chem_group must be an int >= 1 "
f"(got {blocks_per_chem_group})")
max_ext = seed_size - 1 # -1 => start seed already has size 1
ref_chem_groups = copy.deepcopy(the_greed.ref_chem_groups)
mdl_chem_groups = copy.deepcopy(the_greed.mdl_chem_groups)
mapping = dict()
something_happened = True
while something_happened:
something_happened = False
starting_blocks = list()
for ref_chains, mdl_chains in zip(ref_chem_groups, mdl_chem_groups):
if len(mdl_chains) == 0:
continue # nothing to map
ref_chains_copy = list(ref_chains)
for i in range(blocks_per_chem_group):
if len(ref_chains_copy) == 0:
break
seeds = list()
for ref_ch in ref_chains_copy:
seeds += [(ref_ch, mdl_ch) for mdl_ch in mdl_chains]
# extend starting seeds to *seed_size* and retain best scoring
# block for further extension
best_score = -1.0
best_mapping = None
best_seed = None
for s in seeds:
seed = dict(mapping)
seed.update({s[0]: s[1]})
seed = the_greed.ExtendMapping(seed, max_ext = max_ext)
seed_lddt = the_greed.lDDTFromFlatMap(seed)
if seed_lddt > best_score:
best_score = seed_lddt
best_mapping = seed
best_seed = s
if best_mapping != None:
starting_blocks.append(best_mapping)
if best_seed[0] in ref_chains_copy:
# remove that ref chain to enforce diversity
ref_chains_copy.remove(best_seed[0])
# fully expand initial starting blocks
best_lddt = 0.0
best_mapping = None
for seed in starting_blocks:
seed = the_greed.ExtendMapping(seed)
seed_lddt = the_greed.lDDTFromFlatMap(seed)
if seed_lddt > best_lddt:
best_lddt = seed_lddt
best_mapping = seed
if best_lddt == 0.0:
break # no proper mapping found anymore
something_happened = True
mapping.update(best_mapping)
for ref_ch, mdl_ch in best_mapping.items():
for group_idx in range(len(ref_chem_groups)):
if ref_ch in ref_chem_groups[group_idx]:
ref_chem_groups[group_idx].remove(ref_ch)
if mdl_ch in mdl_chem_groups[group_idx]:
mdl_chem_groups[group_idx].remove(mdl_ch)
# translate mapping format and return
final_mapping = list()
for ref_chains in the_greed.ref_chem_groups:
mapped_mdl_chains = list()
for ref_ch in ref_chains:
if ref_ch in mapping:
mapped_mdl_chains.append(mapping[ref_ch])
else:
mapped_mdl_chains.append(None)
final_mapping.append(mapped_mdl_chains)
return final_mapping
class _QSScoreGreedySearcher(qsscore.QSScorer):
def __init__(self, ref, mdl, ref_chem_groups, mdl_chem_groups,
ref_mdl_alns, contact_d = 12.0,
steep_opt_rate = None, greedy_prune_contact_map=False):
""" Greedy extension of already existing but incomplete chain mappings
"""
super().__init__(ref, ref_chem_groups, mdl, ref_mdl_alns,
contact_d = contact_d)
self.ref = ref
self.mdl = mdl
self.ref_mdl_alns = ref_mdl_alns
self.steep_opt_rate = steep_opt_rate
if greedy_prune_contact_map:
self.neighbors = {k: set() for k in self.qsent1.chain_names}
for p in self.qsent1.interacting_chains:
if np.count_nonzero(self.qsent1.PairDist(p[0], p[1])<=8) >= 3:
self.neighbors[p[0]].add(p[1])
self.neighbors[p[1]].add(p[0])
self.mdl_neighbors = {k: set() for k in self.qsent2.chain_names}
for p in self.qsent2.interacting_chains:
if np.count_nonzero(self.qsent2.PairDist(p[0], p[1])<=8) >= 3:
self.mdl_neighbors[p[0]].add(p[1])
self.mdl_neighbors[p[1]].add(p[0])
else:
self.neighbors = {k: set() for k in self.qsent1.chain_names}
for p in self.qsent1.interacting_chains:
self.neighbors[p[0]].add(p[1])
self.neighbors[p[1]].add(p[0])
self.mdl_neighbors = {k: set() for k in self.qsent2.chain_names}
for p in self.qsent2.interacting_chains:
self.mdl_neighbors[p[0]].add(p[1])
self.mdl_neighbors[p[1]].add(p[0])
assert(len(ref_chem_groups) == len(mdl_chem_groups))
self.ref_chem_groups = ref_chem_groups
self.mdl_chem_groups = mdl_chem_groups
self.ref_ch_group_mapper = dict()
self.mdl_ch_group_mapper = dict()
for g_idx, (ref_g, mdl_g) in enumerate(zip(ref_chem_groups,
mdl_chem_groups)):
for ch in ref_g:
self.ref_ch_group_mapper[ch] = g_idx
for ch in mdl_g:
self.mdl_ch_group_mapper[ch] = g_idx
# cache for lDDT based single chain conserved contacts
# used to identify starting points for further extension by QS score
# key: tuple (ref_ch, mdl_ch) value: number of conserved contacts
self.single_chain_scorer = dict()
self.single_chain_cache = dict()
for ch in self.ref.chains:
single_chain_ref = _CSel(self.ref, [ch.GetName()])
self.single_chain_scorer[ch.GetName()] = \
lddt.lDDTScorer(single_chain_ref, bb_only = True)
def SCCounts(self, ref_ch, mdl_ch):
if not (ref_ch, mdl_ch) in self.single_chain_cache:
alns = dict()
alns[mdl_ch] = self.ref_mdl_alns[(ref_ch, mdl_ch)]
mdl_sel = _CSel(self.mdl, [mdl_ch])
s = self.single_chain_scorer[ref_ch]
_,_,_,conserved,_,_,_ = s.lDDT(mdl_sel,
residue_mapping=alns,
return_dist_test=True,
no_interchain=True,
chain_mapping={mdl_ch: ref_ch},
check_resnames=False)
self.single_chain_cache[(ref_ch, mdl_ch)] = conserved
return self.single_chain_cache[(ref_ch, mdl_ch)]
def ExtendMapping(self, mapping, max_ext = None):
if len(mapping) == 0:
raise RuntimError("Mapping must contain a starting point")
# Ref chains onto which we can map. The algorithm starts with a mapping
# on ref_ch. From there we can start to expand to connected neighbors.
# All neighbors that we can reach from the already mapped chains are
# stored in this set which will be updated during runtime
map_targets = set()
for ref_ch in mapping.keys():
map_targets.update(self.neighbors[ref_ch])
# remove the already mapped chains
for ref_ch in mapping.keys():
map_targets.discard(ref_ch)
if len(map_targets) == 0:
return mapping # nothing to extend
# keep track of what model chains are not yet mapped for each chem group
free_mdl_chains = list()
for chem_group in self.mdl_chem_groups:
tmp = [x for x in chem_group if x not in mapping.values()]
free_mdl_chains.append(set(tmp))
# keep track of what ref chains got a mapping
newly_mapped_ref_chains = list()
something_happened = True
while something_happened:
something_happened=False
if self.steep_opt_rate is not None:
n_chains = len(newly_mapped_ref_chains)
if n_chains > 0 and n_chains % self.steep_opt_rate == 0:
mapping = self._SteepOpt(mapping, newly_mapped_ref_chains)
if max_ext is not None and len(newly_mapped_ref_chains) >= max_ext:
break
score_result = self.FromFlatMapping(mapping)
old_score = score_result.QS_global
nominator = score_result.weighted_scores
denominator = score_result.weight_sum + score_result.weight_extra_all
max_diff = 0.0
max_mapping = None
for ref_ch in map_targets:
chem_group_idx = self.ref_ch_group_mapper[ref_ch]
for mdl_ch in free_mdl_chains[chem_group_idx]:
# we're not computing full QS-score here, we directly hack
# into the QS-score formula to compute a diff
nominator_diff = 0.0
denominator_diff = 0.0
for neighbor in self.neighbors[ref_ch]:
if neighbor in mapping and mapping[neighbor] in \
self.mdl_neighbors[mdl_ch]:
# it's a newly added interface if (ref_ch, mdl_ch)
# are added to mapping
int1 = (ref_ch, neighbor)
int2 = (mdl_ch, mapping[neighbor])
a, b, c, d = self._MappedInterfaceScores(int1, int2)
nominator_diff += a # weighted_scores
denominator_diff += b # weight_sum
denominator_diff += d # weight_extra_all
# the respective interface penalties are subtracted
# from denominator
denominator_diff -= self._InterfacePenalty1(int1)
denominator_diff -= self._InterfacePenalty2(int2)
if nominator_diff > 0:
# Only accept a new solution if its actually connected
# i.e. nominator_diff > 0.
new_nominator = nominator + nominator_diff
new_denominator = denominator + denominator_diff
new_score = 0.0
if new_denominator != 0.0:
new_score = new_nominator/new_denominator
diff = new_score - old_score
if diff > max_diff:
max_diff = diff
max_mapping = (ref_ch, mdl_ch)
if max_mapping is not None:
something_happened = True
# assign new found mapping
mapping[max_mapping[0]] = max_mapping[1]
# add all neighboring chains to map targets as they are now
# reachable
for neighbor in self.neighbors[max_mapping[0]]:
if neighbor not in mapping:
map_targets.add(neighbor)
# remove the ref chain from map targets
map_targets.remove(max_mapping[0])
# remove the mdl chain from free_mdl_chains - its taken...
chem_group_idx = self.ref_ch_group_mapper[max_mapping[0]]
free_mdl_chains[chem_group_idx].remove(max_mapping[1])
# keep track of what ref chains got a mapping
newly_mapped_ref_chains.append(max_mapping[0])
return mapping
def _SteepOpt(self, mapping, chains_to_optimize=None):
# just optimize ALL ref chains if nothing specified
if chains_to_optimize is None:
chains_to_optimize = mapping.keys()
# make sure that we only have ref chains which are actually mapped
ref_chains = [x for x in chains_to_optimize if mapping[x] is not None]
# group ref chains to be optimized into chem groups
tmp = dict()
for ch in ref_chains:
chem_group_idx = self.ref_ch_group_mapper[ch]
if chem_group_idx in tmp:
tmp[chem_group_idx].append(ch)
else:
tmp[chem_group_idx] = [ch]
chem_groups = list(tmp.values())
# try all possible mapping swaps. Swaps that improve the score are
# immediately accepted and we start all over again
score_result = self.FromFlatMapping(mapping)
current_score = score_result.QS_global
something_happened = True
while something_happened:
something_happened = False
for chem_group in chem_groups:
if something_happened:
break
for ch1, ch2 in itertools.combinations(chem_group, 2):
swapped_mapping = dict(mapping)
swapped_mapping[ch1] = mapping[ch2]
swapped_mapping[ch2] = mapping[ch1]
score_result = self.FromFlatMapping(swapped_mapping)
if score_result.QS_global > current_score:
something_happened = True
mapping = swapped_mapping
current_score = score_result.QS_global
break
return mapping
def _QSScoreNaive(trg, mdl, chem_groups, chem_mapping, ref_mdl_alns, contact_d,
n_max_naive):
best_mapping = None
best_score = -1.0
# qs_scorer implements caching, score calculation is thus as fast as it gets
# you'll just hit a wall when the number of possible mappings becomes large
qs_scorer = qsscore.QSScorer(trg, chem_groups, mdl, ref_mdl_alns)
for mapping in _ChainMappings(chem_groups, chem_mapping, n_max_naive):
score_result = qs_scorer.Score(mapping, check=False)
if score_result.QS_global > best_score:
best_mapping = mapping
best_score = score_result.QS_global
return (best_mapping, best_score)
def _QSScoreGreedyFast(the_greed):
something_happened = True
mapping = dict()
while something_happened:
something_happened = False
# search for best scoring starting point, we're using lDDT here
n_best = 0
best_seed = None
seeds = _GetSeeds(the_greed.ref_chem_groups,
the_greed.mdl_chem_groups,
mapped_ref_chains = set(mapping.keys()),
mapped_mdl_chains = set(mapping.values()))
for seed in seeds:
n = the_greed.SCCounts(seed[0], seed[1])
if n > n_best:
n_best = n
best_seed = seed
if n_best == 0:
break # no proper seed found anymore...
# add seed to mapping and start the greed
mapping[best_seed[0]] = best_seed[1]
mapping = the_greed.ExtendMapping(mapping)
something_happened = True
# translate mapping format and return
final_mapping = list()
for ref_chains in the_greed.ref_chem_groups:
mapped_mdl_chains = list()
for ref_ch in ref_chains:
if ref_ch in mapping:
mapped_mdl_chains.append(mapping[ref_ch])
else:
mapped_mdl_chains.append(None)
final_mapping.append(mapped_mdl_chains)
return final_mapping
def _QSScoreGreedyFull(the_greed):
""" Uses each reference chain as starting point for expansion
"""
seeds = _GetSeeds(the_greed.ref_chem_groups, the_greed.mdl_chem_groups)
best_overall_score = -1.0
best_overall_mapping = dict()
for seed in seeds:
# do initial extension
mapping = the_greed.ExtendMapping({seed[0]: seed[1]})
# repeat the process until we have a full mapping
something_happened = True
while something_happened:
something_happened = False
remnant_seeds = _GetSeeds(the_greed.ref_chem_groups,
the_greed.mdl_chem_groups,
mapped_ref_chains = set(mapping.keys()),
mapped_mdl_chains = set(mapping.values()))
if len(remnant_seeds) > 0:
# still more mapping to be done
best_score = -1.0
best_mapping = None
for remnant_seed in remnant_seeds:
tmp_mapping = dict(mapping)
tmp_mapping[remnant_seed[0]] = remnant_seed[1]
tmp_mapping = the_greed.ExtendMapping(tmp_mapping)
score_result = the_greed.FromFlatMapping(tmp_mapping)
if score_result.QS_global > best_score:
best_score = score_result.QS_global
best_mapping = tmp_mapping
if best_mapping is not None:
something_happened = True
mapping = best_mapping
score_result = the_greed.FromFlatMapping(mapping)
if score_result.QS_global > best_overall_score:
best_overall_score = score_result.QS_global
best_overall_mapping = mapping
mapping = best_overall_mapping
# translate mapping format and return
final_mapping = list()
for ref_chains in the_greed.ref_chem_groups:
mapped_mdl_chains = list()
for ref_ch in ref_chains:
if ref_ch in mapping:
mapped_mdl_chains.append(mapping[ref_ch])
else:
mapped_mdl_chains.append(None)
final_mapping.append(mapped_mdl_chains)
return final_mapping
def _QSScoreGreedyBlock(the_greed, seed_size, blocks_per_chem_group):
""" try multiple seeds, i.e. try all ref/mdl chain combinations within the
respective chem groups and compute single chain lDDTs. The
*blocks_per_chem_group* best scoring ones are extend by *seed_size* chains
and the best scoring one with respect to QS score is exhaustively extended.
"""
if seed_size is None or seed_size < 1:
raise RuntimeError(f"seed_size must be an int >= 1 (got {seed_size})")
if blocks_per_chem_group is None or blocks_per_chem_group < 1:
raise RuntimeError(f"blocks_per_chem_group must be an int >= 1 "
f"(got {blocks_per_chem_group})")
max_ext = seed_size - 1 # -1 => start seed already has size 1
ref_chem_groups = copy.deepcopy(the_greed.ref_chem_groups)
mdl_chem_groups = copy.deepcopy(the_greed.mdl_chem_groups)
mapping = dict()
something_happened = True
while something_happened:
something_happened = False
starting_blocks = list()
for ref_chains, mdl_chains in zip(ref_chem_groups, mdl_chem_groups):
if len(mdl_chains) == 0:
continue # nothing to map
ref_chains_copy = list(ref_chains)
for i in range(blocks_per_chem_group):
if len(ref_chains_copy) == 0:
break
seeds = list()
for ref_ch in ref_chains_copy:
seeds += [(ref_ch, mdl_ch) for mdl_ch in mdl_chains]
# extend starting seeds to *seed_size* and retain best scoring block
# for further extension
best_score = -1.0
best_mapping = None
best_seed = None
for s in seeds:
seed = dict(mapping)
seed.update({s[0]: s[1]})
seed = the_greed.ExtendMapping(seed, max_ext = max_ext)
score_result = the_greed.FromFlatMapping(seed)
if score_result.QS_global > best_score:
best_score = score_result.QS_global
best_mapping = seed
best_seed = s
if best_mapping != None:
starting_blocks.append(best_mapping)
if best_seed[0] in ref_chains_copy:
# remove selected ref chain to enforce diversity
ref_chains_copy.remove(best_seed[0])
# fully expand initial starting blocks
best_score = -1.0
best_mapping = None
for seed in starting_blocks:
seed = the_greed.ExtendMapping(seed)
score_result = the_greed.FromFlatMapping(seed)
if score_result.QS_global > best_score:
best_score = score_result.QS_global
best_mapping = seed
if best_mapping is not None and len(best_mapping) > len(mapping):
# this even accepts extensions that lead to no increase in QS-score
# at least they make sense from an lDDT perspective
something_happened = True
mapping.update(best_mapping)
for ref_ch, mdl_ch in best_mapping.items():
for group_idx in range(len(ref_chem_groups)):
if ref_ch in ref_chem_groups[group_idx]:
ref_chem_groups[group_idx].remove(ref_ch)
if mdl_ch in mdl_chem_groups[group_idx]:
mdl_chem_groups[group_idx].remove(mdl_ch)
# translate mapping format and return
final_mapping = list()
for ref_chains in the_greed.ref_chem_groups:
mapped_mdl_chains = list()
for ref_ch in ref_chains:
if ref_ch in mapping:
mapped_mdl_chains.append(mapping[ref_ch])
else:
mapped_mdl_chains.append(None)
final_mapping.append(mapped_mdl_chains)
return final_mapping
def _SingleRigidGDTTS(initial_transforms, initial_mappings, chem_groups,
chem_mapping, trg_group_pos, mdl_group_pos,
single_chain_gdtts_thresh, iterative_superposition,
first_complete, n_trg_chains, n_mdl_chains):
""" Takes initial transforms and sequentially adds chain pairs with
best scoring gdtts that fulfill single_chain_gdtts_thresh. The mapping
from the transform that leads to best overall gdtts score is returned.
Optionally, the first complete mapping, i.e. a mapping that covers all
target chains or all model chains, is returned.
"""
best_mapping = dict()
best_gdt = 0
for transform in initial_transforms:
mapping = dict()
mapped_mdl_chains = set()
gdt = 0.0
for trg_chains, mdl_chains, trg_pos, mdl_pos, in zip(chem_groups,
chem_mapping,
trg_group_pos,
mdl_group_pos):
if len(trg_pos) == 0 or len(mdl_pos) == 0:
continue # cannot compute valid gdt
gdt_scores = list()
t_mdl_pos = list()
for m_pos in mdl_pos:
t_m_pos = geom.Vec3List(m_pos)
t_m_pos.ApplyTransform(transform)
t_mdl_pos.append(t_m_pos)
for t_pos, t in zip(trg_pos, trg_chains):
for t_m_pos, m in zip(t_mdl_pos, mdl_chains):
gdt = t_pos.GetGDTTS(t_m_pos)
if gdt >= single_chain_gdtts_thresh:
gdt_scores.append((gdt, (t,m)))
n_gdt_contacts = 4 * len(trg_pos[0])
gdt_scores.sort(reverse=True)
for item in gdt_scores:
p = item[1]
if p[0] not in mapping and p[1] not in mapped_mdl_chains:
mapping[p[0]] = p[1]
mapped_mdl_chains.add(p[1])
gdt += (item[0] * n_gdt_contacts)
if gdt > best_gdt:
best_gdt = gdt
best_mapping = mapping
if first_complete:
n = len(mapping)
if n == n_mdl_chains or n == n_trg_chains:
break
return best_mapping
def _IterativeRigidGDTTS(initial_transforms, initial_mappings, chem_groups,
chem_mapping, trg_group_pos, mdl_group_pos,
single_chain_gdtts_thresh, iterative_superposition,
first_complete, n_trg_chains, n_mdl_chains):
""" Takes initial transforms and sequentially adds chain pairs with
best scoring gdtts that fulfill single_chain_gdtts_thresh. With each
added chain pair, the transform gets updated. Thus the naming iterative.
The mapping from the initial transform that leads to best overall gdtts
score is returned. Optionally, the first complete mapping, i.e. a mapping
that covers all target chains or all model chains, is returned.
"""
# to directly retrieve positions using chain names
trg_pos_dict = dict()
for trg_pos, trg_chains in zip(trg_group_pos, chem_groups):
for t_pos, t in zip(trg_pos, trg_chains):
trg_pos_dict[t] = t_pos
mdl_pos_dict = dict()
for mdl_pos, mdl_chains in zip(mdl_group_pos, chem_mapping):
for m_pos, m in zip(mdl_pos, mdl_chains):
mdl_pos_dict[m] = m_pos
best_mapping = dict()
best_gdt = 0
for initial_transform, initial_mapping in zip(initial_transforms,
initial_mappings):
mapping = {initial_mapping[0]: initial_mapping[1]}
transform = geom.Mat4(initial_transform)
mapped_trg_pos = geom.Vec3List(trg_pos_dict[initial_mapping[0]])
mapped_mdl_pos = geom.Vec3List(mdl_pos_dict[initial_mapping[1]])
# the following variables contain the chains which are
# available for mapping
trg_chain_groups = [set(group) for group in chem_groups]
mdl_chain_groups = [set(group) for group in chem_mapping]
# search and kick out inital mapping
for group in trg_chain_groups:
if initial_mapping[0] in group:
group.remove(initial_mapping[0])
break
for group in mdl_chain_groups:
if initial_mapping[1] in group:
group.remove(initial_mapping[1])
break
something_happened = True
while something_happened:
# search for best mapping given current transform
something_happened=False
best_sc_mapping = None
best_sc_group_idx = None
best_sc_gdt = 0.0
group_idx = 0
for trg_chains, mdl_chains in zip(trg_chain_groups, mdl_chain_groups):
for t in trg_chains:
t_pos = trg_pos_dict[t]
for m in mdl_chains:
m_pos = mdl_pos_dict[m]
t_m_pos = geom.Vec3List(m_pos)
t_m_pos.ApplyTransform(transform)
gdt = t_pos.GetGDTTS(t_m_pos)
if gdt > single_chain_gdtts_thresh and gdt > best_sc_gdt:
best_sc_gdt = gdt
best_sc_mapping = (t,m)
best_sc_group_idx = group_idx
group_idx += 1
if best_sc_mapping is not None:
something_happened = True
mapping[best_sc_mapping[0]] = best_sc_mapping[1]
mapped_trg_pos.extend(trg_pos_dict[best_sc_mapping[0]])
mapped_mdl_pos.extend(mdl_pos_dict[best_sc_mapping[1]])
trg_chain_groups[best_sc_group_idx].remove(best_sc_mapping[0])
mdl_chain_groups[best_sc_group_idx].remove(best_sc_mapping[1])
transform = _GetTransform(mapped_mdl_pos, mapped_trg_pos,
iterative_superposition)
# compute overall gdt for current transform (non-normalized gdt!!!)
mapped_mdl_pos.ApplyTransform(transform)
gdt = mapped_trg_pos.GetGDTTS(mapped_mdl_pos, norm=False)
if gdt > best_gdt:
best_gdt = gdt
best_mapping = mapping
if first_complete:
n = len(mapping)
if n == n_mdl_chains or n == n_trg_chains:
break
return best_mapping
def _SingleRigidRMSD(initial_transforms, initial_mappings, chem_groups,
chem_mapping, trg_group_pos, mdl_group_pos,
iterative_superposition):
"""
Takes initial transforms and sequentially adds chain pairs with lowest RMSD.
The mapping from the transform that leads to lowest overall RMSD is
returned.
"""
best_mapping = dict()
best_ssd = float("inf") # we're actually going for summed squared distances
# Since all positions have same lengths and we do a
# full mapping, lowest SSD has a guarantee of also
# being lowest RMSD
for transform in initial_transforms:
mapping = dict()
mapped_mdl_chains = set()
ssd = 0.0
for trg_chains, mdl_chains, trg_pos, mdl_pos, in zip(chem_groups,
chem_mapping,
trg_group_pos,
mdl_group_pos):
if len(trg_pos) == 0 or len(mdl_pos) == 0:
continue # cannot compute valid rmsd
ssds = list()
t_mdl_pos = list()
for m_pos in mdl_pos:
t_m_pos = geom.Vec3List(m_pos)
t_m_pos.ApplyTransform(transform)
t_mdl_pos.append(t_m_pos)
for t_pos, t in zip(trg_pos, trg_chains):
for t_m_pos, m in zip(t_mdl_pos, mdl_chains):
ssd = t_pos.GetSummedSquaredDistances(t_m_pos)
ssds.append((ssd, (t,m)))
ssds.sort()
for item in ssds:
p = item[1]
if p[0] not in mapping and p[1] not in mapped_mdl_chains:
mapping[p[0]] = p[1]
mapped_mdl_chains.add(p[1])
ssd += item[0]
if ssd < best_ssd:
best_ssd = ssd
best_mapping = mapping
return best_mapping
def _IterativeRigidRMSD(initial_transforms, initial_mappings, chem_groups,
chem_mapping, trg_group_pos, mdl_group_pos,
iterative_superposition):
""" Takes initial transforms and sequentially adds chain pairs with
lowest RMSD. With each added chain pair, the transform gets updated.
Thus the naming iterative. The mapping from the initial transform that
leads to best overall RMSD score is returned.
"""
# to directly retrieve positions using chain names
trg_pos_dict = dict()
for trg_pos, trg_chains in zip(trg_group_pos, chem_groups):
for t_pos, t in zip(trg_pos, trg_chains):
trg_pos_dict[t] = t_pos
mdl_pos_dict = dict()
for mdl_pos, mdl_chains in zip(mdl_group_pos, chem_mapping):
for m_pos, m in zip(mdl_pos, mdl_chains):
mdl_pos_dict[m] = m_pos
best_mapping = dict()
best_rmsd = float("inf")
for initial_transform, initial_mapping in zip(initial_transforms,
initial_mappings):
mapping = {initial_mapping[0]: initial_mapping[1]}
transform = geom.Mat4(initial_transform)
mapped_trg_pos = geom.Vec3List(trg_pos_dict[initial_mapping[0]])
mapped_mdl_pos = geom.Vec3List(mdl_pos_dict[initial_mapping[1]])
# the following variables contain the chains which are
# available for mapping
trg_chain_groups = [set(group) for group in chem_groups]
mdl_chain_groups = [set(group) for group in chem_mapping]
# search and kick out inital mapping
for group in trg_chain_groups:
if initial_mapping[0] in group:
group.remove(initial_mapping[0])
break
for group in mdl_chain_groups:
if initial_mapping[1] in group:
group.remove(initial_mapping[1])
break
something_happened = True
while something_happened:
# search for best mapping given current transform
something_happened=False
best_sc_mapping = None
best_sc_group_idx = None
best_sc_rmsd = float("inf")
group_idx = 0
for trg_chains, mdl_chains in zip(trg_chain_groups, mdl_chain_groups):
for t in trg_chains:
t_pos = trg_pos_dict[t]
for m in mdl_chains:
m_pos = mdl_pos_dict[m]
t_m_pos = geom.Vec3List(m_pos)
t_m_pos.ApplyTransform(transform)
rmsd = t_pos.GetRMSD(t_m_pos)
if rmsd < best_sc_rmsd:
best_sc_rmsd = rmsd
best_sc_mapping = (t,m)
best_sc_group_idx = group_idx
group_idx += 1
if best_sc_mapping is not None:
something_happened = True
mapping[best_sc_mapping[0]] = best_sc_mapping[1]
mapped_trg_pos.extend(trg_pos_dict[best_sc_mapping[0]])
mapped_mdl_pos.extend(mdl_pos_dict[best_sc_mapping[1]])
trg_chain_groups[best_sc_group_idx].remove(best_sc_mapping[0])
mdl_chain_groups[best_sc_group_idx].remove(best_sc_mapping[1])
transform = _GetTransform(mapped_mdl_pos, mapped_trg_pos,
iterative_superposition)
# compute overall RMSD for current transform
mapped_mdl_pos.ApplyTransform(transform)
rmsd = mapped_trg_pos.GetRMSD(mapped_mdl_pos)
if rmsd < best_rmsd:
best_rmsd = rmsd
best_mapping = mapping
return best_mapping
def _GetRefPos(trg, mdl, trg_msas, mdl_alns, max_pos = None):
""" Extracts reference positions which are present in trg and mdl
"""
# select only backbone atoms, makes processing simpler later on
# (just select res.atoms[0].GetPos() as ref pos)
bb_trg = trg.Select("aname=\"CA\",\"C3'\"")
bb_mdl = mdl.Select("aname=\"CA\",\"C3'\"")
# mdl_alns are pairwise, let's construct MSAs
mdl_msas = list()
for aln_list in mdl_alns:
if len(aln_list) > 0:
tmp = aln_list[0].GetSequence(0)
ref_seq = seq.CreateSequence(tmp.GetName(), tmp.GetGaplessString())
mdl_msas.append(seq.alg.MergePairwiseAlignments(aln_list, ref_seq))
else:
mdl_msas.append(seq.CreateAlignment())
trg_pos = list()
mdl_pos = list()
for trg_msa, mdl_msa in zip(trg_msas, mdl_msas):
if mdl_msa.GetCount() > 0:
# make sure they have the same ref sequence (should be a given...)
assert(trg_msa.GetSequence(0).GetGaplessString() == \
mdl_msa.GetSequence(0).GetGaplessString())
else:
# if mdl_msa is empty, i.e. no model chain maps to the chem group
# represented by trg_msa, we just continue. The result will be
# empty position lists added to trg_pos and mdl_pos.
pass
# check which columns in MSAs are fully covered (indices relative to
# first sequence)
trg_indices = _GetFullyCoveredIndices(trg_msa)
mdl_indices = _GetFullyCoveredIndices(mdl_msa)
# get indices where both, mdl and trg, are fully covered
indices = sorted(list(trg_indices.intersection(mdl_indices)))
# subsample if necessary
if max_pos is not None and len(indices) > max_pos:
step = int(len(indices)/max_pos)
indices = [indices[i] for i in range(0, len(indices), step)]
# translate to column indices in the respective MSAs
trg_indices = _RefIndicesToColumnIndices(trg_msa, indices)
mdl_indices = _RefIndicesToColumnIndices(mdl_msa, indices)
# extract positions
trg_pos.append(list())
mdl_pos.append(list())
for s_idx in range(trg_msa.GetCount()):
trg_pos[-1].append(_ExtractMSAPos(trg_msa, s_idx, trg_indices,
bb_trg))
# first seq in mdl_msa is ref sequence in trg and does not belong to mdl
for s_idx in range(1, mdl_msa.GetCount()):
mdl_pos[-1].append(_ExtractMSAPos(mdl_msa, s_idx, mdl_indices,
bb_mdl))
return (trg_pos, mdl_pos)
def _GetFullyCoveredIndices(msa):
""" Helper for _GetRefPos
Returns a set containing the indices relative to first sequence in msa which
are fully covered in all other sequences
--AA-A-A
-BBBB-BB
CCCC-C-C
=> (0,1,3)
"""
indices = set()
ref_idx = 0
for col in msa:
if sum([1 for olc in col if olc != '-']) == col.GetRowCount():
indices.add(ref_idx)
if col[0] != '-':
ref_idx += 1
return indices
def _RefIndicesToColumnIndices(msa, indices):
""" Helper for _GetRefPos
Returns a list of mapped indices. indices refer to non-gap one letter
codes in the first msa sequence. The returnes mapped indices are translated
to the according msa column indices
"""
ref_idx = 0
mapping = dict()
for col_idx, col in enumerate(msa):
if col[0] != '-':
mapping[ref_idx] = col_idx
ref_idx += 1
return [mapping[i] for i in indices]
def _ExtractMSAPos(msa, s_idx, indices, view):
""" Helper for _GetRefPos
Returns a geom.Vec3List containing positions refering to given msa sequence.
=> Chain with corresponding name is mapped onto sequence and the position of
the first atom of each residue specified in indices is extracted.
Indices refers to column indices in msa!
"""
s = msa.GetSequence(s_idx)
s_v = _CSel(view, [s.GetName()])
# sanity check
assert(len(s.GetGaplessString()) == len(s_v.residues))
residue_idx = [s.GetResidueIndex(i) for i in indices]
return geom.Vec3List([s_v.residues[i].atoms[0].pos for i in residue_idx])
def _NChemGroupMappings(ref_chains, mdl_chains):
""" Number of mappings within one chem group
:param ref_chains: Reference chains
:type ref_chains: :class:`list` of :class:`str`
:param mdl_chains: Model chains that are mapped onto *ref_chains*
:type mdl_chains: :class:`list` of :class:`str`
:returns: Number of possible mappings of *mdl_chains* onto *ref_chains*
"""
n_ref = len(ref_chains)
n_mdl = len(mdl_chains)
if n_ref == n_mdl:
return factorial(n_ref)
elif n_ref > n_mdl:
n_choose_k = binom(n_ref, n_mdl)
return n_choose_k * factorial(n_mdl)
else:
n_choose_k = binom(n_mdl, n_ref)
return n_choose_k * factorial(n_ref)
def _NMappings(ref_chains, mdl_chains):
""" Number of mappings for a full chem mapping
:param ref_chains: Chem groups of reference
:type ref_chains: :class:`list` of :class:`list` of :class:`str`
:param mdl_chains: Model chains that map onto those chem groups
:type mdl_chains: :class:`list` of :class:`list` of :class:`str`
:returns: Number of possible mappings of *mdl_chains* onto *ref_chains*
"""
assert(len(ref_chains) == len(mdl_chains))
n = 1
for a,b in zip(ref_chains, mdl_chains):
n *= _NChemGroupMappings(a,b)
return n
def _NMappingsWithin(ref_chains, mdl_chains, max_mappings):
""" Check whether total number of mappings is smaller than given maximum
In principle the same as :func:`_NMappings` but it stops as soon as the
maximum is hit.
:param ref_chains: Chem groups of reference
:type ref_chains: :class:`list` of :class:`list` of :class:`str`
:param mdl_chains: Model chains that map onto those chem groups
:type mdl_chains: :class:`list` of :class:`list` of :class:`str`
:param max_mappings: Number of max allowed mappings
:returns: Whether number of possible mappings of *mdl_chains* onto
*ref_chains* is below or equal *max_mappings*.
"""
assert(len(ref_chains) == len(mdl_chains))
n = 1
for a,b in zip(ref_chains, mdl_chains):
n *= _NChemGroupMappings(a,b)
if n > max_mappings:
return False
return True
def _RefSmallerGenerator(ref_chains, mdl_chains):
""" Returns all possible ways to map mdl_chains onto ref_chains
Specific for the case where len(ref_chains) < len(mdl_chains)
"""
for c in itertools.combinations(mdl_chains, len(ref_chains)):
for p in itertools.permutations(c):
yield list(p)
def _RefLargerGenerator(ref_chains, mdl_chains):
""" Returns all possible ways to map mdl_chains onto ref_chains
Specific for the case where len(ref_chains) > len(mdl_chains)
Ref chains without mapped mdl chain are assigned None
"""
n_ref = len(ref_chains)
n_mdl = len(mdl_chains)
for c in itertools.combinations(range(n_ref), n_mdl):
for p in itertools.permutations(mdl_chains):
ret_list = [None] * n_ref
for idx, ch in zip(c, p):
ret_list[idx] = ch
yield ret_list
def _RefEqualGenerator(ref_chains, mdl_chains):
""" Returns all possible ways to map mdl_chains onto ref_chains
Specific for the case where len(ref_chains) == len(mdl_chains)
"""
for p in itertools.permutations(mdl_chains):
yield list(p)
def _ConcatIterators(iterators):
for item in itertools.product(*iterators):
yield list(item)
def _ChainMappings(ref_chains, mdl_chains, n_max=None):
"""Returns all possible ways to map *mdl_chains* onto fixed *ref_chains*
:param ref_chains: List of list of chemically equivalent chains in reference
:type ref_chains: :class:`list` of :class:`list`
:param mdl_chains: Equally long list of list of chemically equivalent chains
in model that map on those ref chains.
:type mdl_chains: :class:`list` of :class:`list`
:param n_max: Aborts and raises :class:`RuntimeError` if max number of
mappings is above this threshold.
:type n_max: :class:`int`
:returns: Iterator over all possible mappings of *mdl_chains* onto fixed
*ref_chains*. Potentially contains None as padding when number of
model chains for a certain mapping is smaller than the according
reference chains.
Example: _ChainMappings([['A', 'B', 'C'], ['D', 'E']],
[['x', 'y'], ['i', 'j']])
gives an iterator over: [[['x', 'y', None], ['i', 'j']],
[['x', 'y', None], ['j', 'i']],
[['y', 'x', None], ['i', 'j']],
[['y', 'x', None], ['j', 'i']],
[['x', None, 'y'], ['i', 'j']],
[['x', None, 'y'], ['j', 'i']],
[['y', None, 'x'], ['i', 'j']],
[['y', None, 'x'], ['j', 'i']],
[[None, 'x', 'y'], ['i', 'j']],
[[None, 'x', 'y'], ['j', 'i']],
[[None, 'y', 'x'], ['i', 'j']],
[[None, 'y', 'x'], ['j', 'i']]]
"""
assert(len(ref_chains) == len(mdl_chains))
if n_max is not None:
if not _NMappingsWithin(ref_chains, mdl_chains, n_max):
raise RuntimeError(f"Too many mappings. Max allowed: {n_max}")
# one iterator per mapping representing all mdl combinations relative to
# reference
iterators = list()
for ref, mdl in zip(ref_chains, mdl_chains):
if len(ref) == 0:
raise RuntimeError("Expext at least one chain in ref chem group")
if len(ref) == len(mdl):
iterators.append(_RefEqualGenerator(ref, mdl))
elif len(ref) < len(mdl):
iterators.append(_RefSmallerGenerator(ref, mdl))
else:
iterators.append(_RefLargerGenerator(ref, mdl))
return _ConcatIterators(iterators)
def _GetTransform(pos_one, pos_two, iterative):
""" Computes minimal RMSD superposition for pos_one onto pos_two
:param pos_one: Positions that should be superposed onto *pos_two*
:type pos_one: :class:`geom.Vec3List`
:param pos_two: Reference positions
:type pos_two: :class:`geom.Vec3List`
:iterative: Whether iterative superposition should be used. Iterative
potentially raises, uses standard superposition as fallback.
:type iterative: :class:`bool`
:returns: Transformation matrix to superpose *pos_one* onto *pos_two*
:rtype: :class:`geom.Mat4`
"""
res = None
if iterative:
try:
res = mol.alg.IterativeSuperposeSVD(pos_one, pos_two)
except:
pass # triggers fallback below
if res is None:
res = mol.alg.SuperposeSVD(pos_one, pos_two)
return res.transformation
# specify public interface
__all__ = ('ChainMapper', 'ReprResult', 'MappingResult')