stage/lib64/python2.7/site-packages/ost/mol/alg/qsscoring.py

Go to the documentation of this file.

"""
Scoring of quaternary structures as in Martino's 2017 paper.

.. note ::

  Requirements for use:

  - A default :class:`compound library <ost.conop.CompoundLib>` must be defined
    and accessible via :func:`~ost.conop.GetDefaultLib`. This is set by default
    when executing scripts with ``ost``. Otherwise, you must set this with
    :func:`~ost.conop.SetDefaultLib`.
  - ClustalW must be installed (unless you provide chain mappings)
  - Python modules `numpy` and `scipy` must be installed and available
    (e.g. use ``pip install scipy numpy``)

Authors: Gerardo Tauriello, Martino Bertoni
"""

from ost import mol, geom, conop, seq, settings
from ost import LogError, LogWarning, LogScript, LogInfo, LogVerbose, LogDebug
from ost.bindings.clustalw import ClustalW
import numpy as np
from scipy.misc import factorial
from scipy.special import binom
from scipy.cluster.hierarchy import fclusterdata
import itertools

###############################################################################
# QS scoring
###############################################################################

class QSscoreError(Exception):
  """Exception to be raised for "acceptable" exceptions in QS scoring.

  Those are cases we might want to capture for default behavior.
  """
  pass

class QSscorer:
  """Object to compute QS scores.

  Simple usage without any precomputed contacts, symmetries and mappings:

  .. code-block:: python

    import ost
    from ost.mol.alg import qsscoring

    # load two biounits to compare
    ent_full = ost.io.LoadPDB('3ia3', remote=True)
    ent_1 = ent_full.Select('cname=A,D')
    ent_2 = ent_full.Select('cname=B,C')
    # get score
    ost.PushVerbosityLevel(3)
    try:
      qs_scorer = qsscoring.QSscorer(ent_1, ent_2)
      ost.LogScript('QSscore:', str(qs_scorer.global_score))
      ost.LogScript('Chain mapping used:', str(qs_scorer.chain_mapping))
    except qsscoring.QSscoreError as ex:
      # default handling: report failure and set score to 0
      ost.LogError('QSscore failed:', str(ex))
      qs_score = 0

  For maximal performance when computing QS scores of the same entity with many
  others, it is advisable to construct and reuse :class:`QSscoreEntity` objects.

  Any known / precomputed information can be filled into the appropriate
  attribute here (no checks done!). Otherwise most quantities are computed on
  first access and cached (lazy evaluation). Setters are provided to set values
  with extra checks (e.g. :func:`SetSymmetries`).

  All necessary seq. alignments are done by global BLOSUM62-based alignment. A
  multiple sequence alignment is performed with ClustalW unless
  :attr:`chain_mapping` is provided manually. You will need to have an
  executable ``clustalw`` or ``clustalw2`` in your ``PATH`` or you must set
  :attr:`clustalw_bin` accordingly. Otherwise an exception
  (:class:`ost.settings.FileNotFound`) is thrown.

  Formulas for QS scores:

  ::
  
    - QS_best = weighted_scores / (weight_sum + weight_extra_mapped)
    - QS_global = weighted_scores / (weight_sum + weight_extra_all)
    -> weighted_scores = sum(w(min(d1,d2)) * (1 - abs(d1-d2)/12)) for shared
    -> weight_sum = sum(w(min(d1,d2))) for shared
    -> weight_extra_mapped = sum(w(d)) for all mapped but non-shared
    -> weight_extra_all = sum(w(d)) for all non-shared
    -> w(d) = 1 if d <= 5, exp(-2 * ((d-5.0)/4.28)^2) else

  :param ent_1: First structure to be scored.
  :type ent_1:  :class:`QSscoreEntity`, :class:`~ost.mol.EntityHandle` or
                :class:`~ost.mol.EntityView`
  :param ent_2: Second structure to be scored.
  :type ent_2:  :class:`QSscoreEntity`, :class:`~ost.mol.EntityHandle` or
                :class:`~ost.mol.EntityView`

  :raises: :class:`QSscoreError` if input structures are invalid or are monomers
           or have issues that make it impossible for a QS score to be computed.

  .. attribute:: qs_ent_1

    :class:`QSscoreEntity` object for *ent_1* given at construction.
    If entity names (:attr:`~QSscoreEntity.original_name`) are not unique, we
    set it to 'pdb_1' using :func:`~QSscoreEntity.SetName`.

  .. attribute:: qs_ent_2

    :class:`QSscoreEntity` object for *ent_2* given at construction.
    If entity names (:attr:`~QSscoreEntity.original_name`) are not unique, we
    set it to 'pdb_2' using :func:`~QSscoreEntity.SetName`.

  .. attribute:: calpha_only

    True if any of the two structures is CA-only (after cleanup).

    :type: :class:`bool`

  .. attribute:: max_ca_per_chain_for_cm

    Maximal number of CA atoms to use in each chain to determine chain mappings.
    Setting this to -1 disables the limit. Limiting it speeds up determination
    of symmetries and chain mappings. By default it is set to 100.

    :type: :class:`int`
  """
  def __init__(self, ent_1, ent_2):
    # generate QSscoreEntity objects?
    if isinstance(ent_1, QSscoreEntity):
      self.qs_ent_1 = ent_1
    else:
      self.qs_ent_1 = QSscoreEntity(ent_1)
    if isinstance(ent_2, QSscoreEntity):
      self.qs_ent_2 = ent_2
    else:
      self.qs_ent_2 = QSscoreEntity(ent_2)
    # check validity of inputs
    if not self.qs_ent_1.is_valid or not self.qs_ent_2.is_valid:
      raise QSscoreError("Invalid input in QSscorer!")
    # set names for structures
    if self.qs_ent_1.original_name == self.qs_ent_2.original_name:
      self.qs_ent_1.SetName('pdb_1')
      self.qs_ent_2.SetName('pdb_2')
    else:
      self.qs_ent_1.SetName(self.qs_ent_1.original_name)
      self.qs_ent_2.SetName(self.qs_ent_2.original_name)
    # set other public attributes
    self.calpha_only = self.qs_ent_1.calpha_only or self.qs_ent_2.calpha_only
    self.max_ca_per_chain_for_cm = 100
    # init cached stuff
    self._chem_mapping = None
    self._ent_to_cm_1 = None
    self._ent_to_cm_2 = None
    self._symm_1 = None
    self._symm_2 = None
    self._chain_mapping = None
    self._alignments = None
    self._mapped_residues = None
    self._global_score = None
    self._best_score = None
    self._superposition = None
    self._lddt_score = None
    self._lddt_mdl = None
    self._lddt_ref = None
    self._clustalw_bin = None

  @property
  def chem_mapping(self):
    """Inter-complex mapping of chemical groups.

    Each group (see :attr:`QSscoreEntity.chem_groups`) is mapped according to
    highest sequence identity. Alignment is between longest sequences in groups.

    Limitations:

    - If different numbers of groups, we map only the groups for the complex
      with less groups (rest considered unmapped and shown as warning)
    - The mapping is forced: the "best" mapping will be chosen independently of
      how low the seq. identity may be

    :getter: Computed on first use (cached)
    :type: :class:`dict` with key = :class:`tuple` of chain names in
           :attr:`qs_ent_1` and value = :class:`tuple` of chain names in
           :attr:`qs_ent_2`.

    :raises: :class:`QSscoreError` if we end up having less than 2 chains for
             either entity in the mapping (can happen if chains do not have CA
             atoms).
    """
    if self._chem_mapping is None:
      self._chem_mapping = _GetChemGroupsMapping(self.qs_ent_1, self.qs_ent_2)
    return self._chem_mapping

  @property
  def ent_to_cm_1(self):
    """Subset of :attr:`qs_ent_1` used to compute chain mapping and symmetries.

    Properties:

    - Includes only residues aligned according to :attr:`chem_mapping`
    - Includes only 1 CA atom per residue
    - Has at least 5 and at most :attr:`max_ca_per_chain_for_cm` atoms per chain
    - All chains of the same chemical group have the same number of atoms
      (also in :attr:`ent_to_cm_2` according to :attr:`chem_mapping`)
    - All chains appearing in :attr:`chem_mapping` appear in this entity
      (so the two can be safely used together)

    This entity might be transformed (i.e. all positions rotated/translated by
    same transformation matrix) if this can speed up computations. So do not
    assume fixed global positions (but relative distances will remain fixed).

    :getter: Computed on first use (cached)
    :type: :class:`~ost.mol.EntityHandle`

    :raises: :class:`QSscoreError` if any chain ends up having less than 5 res.
    """
    if self._ent_to_cm_1 is None:
      self._ComputeAlignedEntities()
    return self._ent_to_cm_1

  @property
  def ent_to_cm_2(self):
    """Subset of :attr:`qs_ent_1` used to compute chain mapping and symmetries
    (see :attr:`ent_to_cm_1` for details).
    """
    if self._ent_to_cm_2 is None:
      self._ComputeAlignedEntities()
    return self._ent_to_cm_2

  @property
  def symm_1(self):
    """Symmetry groups for :attr:`qs_ent_1` used to speed up chain mapping.

    This is a list of chain-lists where each chain-list can be used reconstruct
    the others via cyclic C or dihedral D symmetry. The first chain-list is used
    as a representative symmetry group. For heteromers, the group-members must
    contain all different seqres in oligomer.

    Example: symm. groups [(A,B,C), (D,E,F), (G,H,I)] means that there are
    symmetry transformations to get (D,E,F) and (G,H,I) from (A,B,C).

    Properties:

    - All symmetry group tuples have the same length (num. of chains)
    - All chains in :attr:`ent_to_cm_1` appear (w/o duplicates)
    - For heteros: symmetry group tuples have all different chem. groups
    - Trivial symmetry group = one tuple with all chains (used if inconsistent
      data provided or if no symmetry is found)
    - Either compatible to :attr:`symm_2` or trivial symmetry groups used.
      Compatibility requires same lengths of symmetry group tuples and it must
      be possible to get an overlap (80% of residues covered within 6 A of a
      (chem. mapped) chain) of all chains in representative symmetry groups by
      superposing one pair of chains.

    :getter: Computed on first use (cached)
    :type: :class:`list` of :class:`tuple` of :class:`str` (chain names)
    """
    if self._symm_1 is None:
      self._ComputeSymmetry()
    return self._symm_1

  @property
  def symm_2(self):
    """Symmetry groups for :attr:`qs_ent_2` (see :attr:`symm_1` for details)."""
    if self._symm_2 is None:
      self._ComputeSymmetry()
    return self._symm_2

  def SetSymmetries(self, symm_1, symm_2):
    """Set user-provided symmetry groups.

    These groups are restricted to chain names appearing in :attr:`ent_to_cm_1`
    and :attr:`ent_to_cm_2` respectively. They are only valid if they cover all
    chains and both *symm_1* and *symm_2* have same lengths of symmetry group
    tuples. Otherwise trivial symmetry group used (see :attr:`symm_1`).

    :param symm_1: Value to set for :attr:`symm_1`.
    :param symm_2: Value to set for :attr:`symm_2`.
    """
    # restrict chain names
    self._symm_1 = _CleanUserSymmetry(symm_1, self.ent_to_cm_1)
    self._symm_2 = _CleanUserSymmetry(symm_2, self.ent_to_cm_2)
    # check that we have reasonable symmetries set (fallback: all chains)
    if not _AreValidSymmetries(self._symm_1, self._symm_2):
      self._symm_1 = [tuple(ch.name for ch in self.ent_to_cm_1.chains)]
      self._symm_2 = [tuple(ch.name for ch in self.ent_to_cm_2.chains)]

  @property
  def chain_mapping(self):
    """Mapping from :attr:`ent_to_cm_1` to :attr:`ent_to_cm_2`.

    Properties:

    - Mapping is between chains of same chem. group (see :attr:`chem_mapping`)
    - Each chain can appear only once in mapping
    - All chains of complex with less chains are mapped
    - Symmetry (:attr:`symm_1`, :attr:`symm_2`) is taken into account

    Details on algorithms used to find mapping:

    - We try all pairs of chem. mapped chains within symmetry group and get
      superpose-transformation for them
    - First option: check for "sufficient overlap" of other chain-pairs

      - For each chain-pair defined above: apply superposition to full oligomer
        and map chains based on structural overlap
      - Structural overlap = X% of residues in second oligomer covered within Y
        Angstrom of a (chem. mapped) chain in first oligomer. We successively
        try (X,Y) = (80,4), (40,6) and (20,8) to be less and less strict in
        mapping (warning shown for most permissive one).
      - If multiple possible mappings are found, we choose the one which leads
        to the lowest multi-chain-RMSD given the superposition

    - Fallback option: try all mappings to find minimal multi-chain-RMSD
      (warning shown)

      - For each chain-pair defined above: apply superposition, try all (!)
        possible chain mappings (within symmetry group) and keep mapping with
        lowest multi-chain-RMSD
      - Repeat procedure above to resolve symmetry. Within the symmetry group we
        can use the chain mapping computed before and we just need to find which
        symmetry group in first oligomer maps to which in the second one. We
        again try all possible combinations...
      - Limitations:
        
        - Trying all possible mappings is a combinatorial nightmare (factorial).
          We throw an exception if too many combinations (e.g. octomer vs
          octomer with no usable symmetry)
        - The mapping is forced: the "best" mapping will be chosen independently
          of how badly they fit in terms of multi-chain-RMSD
        - As a result, such a forced mapping can lead to a large range of
          resulting QS scores. An extreme example was observed between 1on3.1
          and 3u9r.1, where :attr:`global_score` can range from 0.12 to 0.43
          for mappings with very similar multi-chain-RMSD.

    :getter: Computed on first use (cached)
    :type: :class:`dict` with key / value = :class:`str` (chain names, key
           for :attr:`ent_to_cm_1`, value for :attr:`ent_to_cm_2`)
    :raises: :class:`QSscoreError` if there are too many combinations to check
             to find a chain mapping.
    """
    if self._chain_mapping is None:
      self._chain_mapping = _GetChainMapping(self.ent_to_cm_1, self.ent_to_cm_2,
                                             self.symm_1, self.symm_2,
                                             self.chem_mapping)
      LogInfo('Mapping found: %s' % str(self._chain_mapping))
    return self._chain_mapping

  @property
  def alignments(self):
    """List of successful sequence alignments using :attr:`chain_mapping`.

    There will be one alignment for each mapped chain and they are ordered by
    their chain names in :attr:`qs_ent_1`.

    The sequences of the alignments have views attached into
    :attr:`QSscoreEntity.ent` of :attr:`qs_ent_1` and :attr:`qs_ent_2`.

    :getter: Computed on first use (cached)
    :type: :class:`list` of :class:`~ost.seq.AlignmentHandle`
    """
    if self._alignments is None:
      self._alignments = _GetMappedAlignments(self.qs_ent_1.ent,
                                              self.qs_ent_2.ent,
                                              self.chain_mapping)
    return self._alignments

  @property
  def mapped_residues(self):
    """Mapping of shared residues in :attr:`alignments`.

    :getter: Computed on first use (cached)
    :type: :class:`dict` *mapped_residues[c1][r1] = r2* with:
           *c1* = Chain name in first entity (= first sequence in aln),
           *r1* = Residue number in first entity,
           *r2* = Residue number in second entity
    """
    if self._mapped_residues is None:
      self._mapped_residues = _GetMappedResidues(self.alignments)
    return self._mapped_residues

  @property
  def global_score(self):
    """QS-score with penalties.
    
    The range of the score is between 0 (i.e. no interface residues are shared
    between biounits) and 1 (i.e. the interfaces are identical).
    
    The global QS-score is computed applying penalties when interface residues
    or entire chains are missing (i.e. anything that is not mapped in
    :attr:`mapped_residues` / :attr:`chain_mapping`) in one of the biounits.

    :getter: Computed on first use (cached)
    :type: :class:`float`
    """
    if self._global_score is None:
      self._ComputeScores()
    return self._global_score

  @property
  def best_score(self):
    """QS-score without penalties.

    Like :attr:`global_score`, but neglecting additional residues or chains in
    one of the biounits (i.e. the score is calculated considering only mapped
    chains and residues).

    :getter: Computed on first use (cached)
    :type: :class:`float`
    """
    if self._best_score is None:
      self._ComputeScores()
    return self._best_score

  @property
  def superposition(self):
    """Superposition result based on shared CA atoms in :attr:`alignments`.

    The superposition can be used to map :attr:`QSscoreEntity.ent` of
    :attr:`qs_ent_1` onto the one of :attr:`qs_ent_2`. Use
    :func:`ost.geom.Invert` if you need the opposite transformation.

    :getter: Computed on first use (cached)
    :type: :class:`ost.mol.alg.SuperpositionResult`
    """
    if self._superposition is None:
      self._superposition = _GetQsSuperposition(self.alignments)
      # report it
      sup_rmsd = self._superposition.rmsd
      cmp_view = self._superposition.view1
      LogInfo('CA RMSD for %s aligned residues on %s chains: %.2f' \
              % (cmp_view.residue_count, cmp_view.chain_count, sup_rmsd))
    return self._superposition

  @property
  def lddt_score(self):
    """The multi-chain lDDT score.

    .. note::

      lDDT is not considering over-prediction (i.e. extra chains) and hence is
      not symmetric. Here, we consider :attr:`qs_ent_1` as the reference and
      :attr:`qs_ent_2` as the model. The alignments from :attr:`alignments` are
      used to map residue numbers and chains.

    The score is computed with OST's :func:`~ost.mol.alg.LocalDistDiffTest`
    function with a single distance threshold of 2 A and an inclusion radius of
    8 A. You can use :attr:`lddt_mdl` and :attr:`lddt_ref` to get entities on
    which you can call any other lDDT function with any other set of parameters.

    :getter: Computed on first use (cached)
    :type: :class:`float`
    """
    if self._lddt_score is None:
      self._ComputeLDDT()
    return self._lddt_score

  @property
  def lddt_mdl(self):
    """The model entity used for lDDT scoring (:attr:`lddt_score`) and annotated
    with local scores.

    Local scores are available as residue properties named 'lddt' and on each
    atom as a B-factor. Only CA atoms are considered if :attr:`calpha_only` is
    True, otherwise this is an all-atom score.
    
    Since, the lDDT computation requires a single chain with mapped residue
    numbering, all chains are appended into a single chain X with unique residue
    numbers according to the column-index in the alignment. The alignments are
    in the same order as they appear in :attr:`alignments`. Additional residues
    are appended at the end of the chain with unique residue numbers.

    :getter: Computed on first use (cached)
    :type: :class:`~ost.mol.EntityHandle`
    """
    if self._lddt_mdl is None:
      self._ComputeLDDT()
    return self._lddt_mdl

  @property
  def lddt_ref(self):
    """The reference entity used for lDDT scoring (:attr:`lddt_score`).

    This is a single chain X with residue numbers matching ones in
    :attr:`lddt_mdl` where aligned and unique numbers for additional residues.

    :getter: Computed on first use (cached)
    :type: :class:`~ost.mol.EntityHandle`
    """
    if self._lddt_ref is None:
      self._ComputeLDDT()
    return self._lddt_ref

  @property
  def clustalw_bin(self):
    """
    Full path to ``clustalw`` or ``clustalw2`` executable to use for multiple
    sequence alignments (unless :attr:`chain_mapping` is provided manually).

    :getter: Located in path on first use (cached)
    :type: :class:`str`
    """
    if self._clustalw_bin is None:
      self._clustalw_bin = settings.Locate(('clustalw', 'clustalw2'))
    return self._clustalw_bin

  ##############################################################################
  # Class internal helpers (anything that doesnt easily work without this class)
  ##############################################################################

  def _ComputeAlignedEntities(self):
    """Fills cached ent_to_cm_1 and ent_to_cm_2."""
    # get aligned residues via MSA
    ev1, ev2 = _GetAlignedResidues(self.qs_ent_1, self.qs_ent_2,
                                   self.chem_mapping,
                                   self.max_ca_per_chain_for_cm,
                                   self.clustalw_bin)
    # extract new entities
    self._ent_to_cm_1 = mol.CreateEntityFromView(ev1, True)
    self._ent_to_cm_2 = mol.CreateEntityFromView(ev2, True)
    # name them
    self._ent_to_cm_1.SetName(self.qs_ent_1.GetName())
    self._ent_to_cm_2.SetName(self.qs_ent_2.GetName())

  def _ComputeSymmetry(self):
    """Fills cached symm_1 and symm_2."""
    # get them
    self._symm_1, self._symm_2 = \
      _FindSymmetry(self.qs_ent_1, self.qs_ent_2, self.ent_to_cm_1,
                    self.ent_to_cm_2, self.chem_mapping)
    # check that we have reasonable symmetries set (fallback: all chains)
    if not _AreValidSymmetries(self._symm_1, self._symm_2):
      self._symm_1 = [tuple(ch.name for ch in self.ent_to_cm_1.chains)]
      self._symm_2 = [tuple(ch.name for ch in self.ent_to_cm_2.chains)]

  def _ComputeScores(self):
    """Fills cached global_score and best_score."""
    # get contacts
    if self.calpha_only:
      contacts_1 = self.qs_ent_1.contacts_ca
      contacts_2 = self.qs_ent_2.contacts_ca
    else:
      contacts_1 = self.qs_ent_1.contacts
      contacts_2 = self.qs_ent_2.contacts
    # get scores
    scores = _GetScores(contacts_1, contacts_2, self.mapped_residues,
                        self.chain_mapping)
    self._best_score = scores[0]
    self._global_score = scores[1]
    # report scores
    LogInfo('QSscore %s, %s: best: %.2f, global: %.2f' \
            % (self.qs_ent_1.GetName(), self.qs_ent_2.GetName(),
               self._best_score, self._global_score))

  def _ComputeLDDT(self):
    """Fills cached lddt_score, lddt_mdl and lddt_ref."""
    LogInfo('Computing lDDT score')
    # check reference and model
    ref, mdl = self.qs_ent_1.ent, self.qs_ent_2.ent
    LogInfo('Reference %s has: %s chains' % (ref.GetName(), ref.chain_count))
    LogInfo('Model %s has: %s chains' % (mdl.GetName(), mdl.chain_count))
    if mdl.chain_count > ref.chain_count:
      LogWarning('MODEL contains more chains than REFERENCE, '
                 'lDDT is not considering them')
    # get single chain reference and model
    self._lddt_ref, self._lddt_mdl = \
      _MergeAlignedChains(self.alignments, ref, mdl, self.calpha_only)
    # score them (mdl and ref changed) and keep results
    self._lddt_score = _ComputeLDDTScore(self._lddt_ref, self._lddt_mdl)


###############################################################################
# Entity with cached entries for QS scoring
###############################################################################

class QSscoreEntity(object):
  """Entity with cached entries for QS scoring.

  Any known / precomputed information can be filled into the appropriate
  attribute here as long as they are labelled as read/write. Otherwise the
  quantities are computed on first access and cached (lazy evaluation). The
  heaviest load is expected when computing :attr:`contacts` and
  :attr:`contacts_ca`.

  :param ent: Entity to be used for QS scoring. A copy of it will be processed.
  :type ent:  :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`

  .. attribute:: is_valid

    True, if successfully initialized. False, if input structure is monomer or
    has less than 2 protein chains with >= 20 residues.

    :type: :class:`bool`

  .. attribute:: original_name

    Name set for *ent* when object was created.

    :type: :class:`str`

  .. attribute:: ent

    Cleaned version of *ent* passed at construction. Hydrogens are removed, the
    entity is processed with a :class:`~ost.conop.RuleBasedProcessor` and chains
    listed in :attr:`removed_chains` have been removed. The name of this entity
    might change during scoring (see :func:`GetName`). Otherwise, this will be
    fixed.

    :type: :class:`~ost.mol.EntityHandle`

  .. attribute:: removed_chains

    Chains removed from *ent* passed at construction. These are ligand and water
    chains as well as small (< 20 res.) peptides or chains with no amino acids
    (determined by chem. type, which is set by rule based processor).

    :type: :class:`list` of :class:`str`

  .. attribute:: calpha_only

    Whether entity is CA-only (i.e. it has 0 CB atoms)

    :type: :class:`bool`
  """
  def __init__(self, ent):
    # copy entity and process/clean it
    self.original_name = ent.GetName()
    ent = mol.CreateEntityFromView(ent.Select('ele!=H and aname!=HN'), True)
    if not conop.GetDefaultLib():
      raise RuntimeError("QSscore computation requires a compound library!")
    pr = conop.RuleBasedProcessor(conop.GetDefaultLib())
    pr.Process(ent)
    self.ent, self.removed_chains, self.calpha_only = _CleanInputEntity(ent)
    # check if it's suitable for QS scoring
    if self.ent.chain_count == 0:
      LogError('Bad input file: ' + ent.GetName() + '. No chains left after '
               'removing water, ligands and small peptides.')
      self.is_valid = False
    elif self.ent.chain_count == 1:
      LogError('Structure ' + ent.GetName() + ' is a monomer. '
               'QSscore is not defined for monomers.')
      self.is_valid = False
    else:
      self.is_valid = True
    # init cached stuff
    self._ca_entity = None
    self._ca_chains = None
    self._alignments = {}
    self._chem_groups = None
    self._angles = {}
    self._axis = {}
    self._contacts = None
    self._contacts_ca = None

  def GetName(self):
    """Wrapper to :func:`~ost.mol.EntityHandle.GetName` of :attr:`ent`.
    This is used to uniquely identify the entity while scoring. The name may
    therefore change while :attr:`original_name` remains fixed.
    """
    # for duck-typing and convenience
    return self.ent.GetName()

  def SetName(self, new_name):
    """Wrapper to :func:`~ost.mol.EntityHandle.SetName` of :attr:`ent`.
    Use this to change unique identifier while scoring (see :func:`GetName`).
    """
    # for duck-typing and convenience
    self.ent.SetName(new_name)

  @property
  def ca_entity(self):
    """
    Reduced representation of :attr:`ent` with only CA atoms.
    This guarantees that each included residue has exactly one atom.

    :getter: Computed on first use (cached)
    :type: :class:`~ost.mol.EntityHandle`
    """
    if self._ca_entity is None:
      self._ca_entity = _GetCAOnlyEntity(self.ent)
    return self._ca_entity

  @property
  def ca_chains(self):
    """
    Map of chain names in :attr:`ent` to sequences with attached view to CA-only
    chains (into :attr:`ca_entity`). Useful for alignments and superpositions.

    :getter: Computed on first use (cached)
    :type: :class:`dict` (key = :class:`str`,
           value = :class:`~ost.seq.SequenceHandle`)
    """
    if self._ca_chains is None:
      self._ca_chains = dict()
      ca_entity = self.ca_entity
      for ch in ca_entity.chains:
        self._ca_chains[ch.name] = seq.SequenceFromChain(ch.name, ch)
    return self._ca_chains

  def GetAlignment(self, c1, c2):
    """Get sequence alignment of chain *c1* with chain *c2*.
    Computed on first use based on :attr:`ca_chains` (cached).

    :param c1: Chain name for first chain to align.
    :type c1:  :class:`str`
    :param c2: Chain name for second chain to align.
    :type c2:  :class:`str`
    :rtype: :class:`~ost.seq.AlignmentHandle` or None if it failed.
    """
    if (c1,c2) not in self._alignments:
      ca_chains = self.ca_chains
      self._alignments[(c1,c2)] = _AlignAtomSeqs(ca_chains[c1], ca_chains[c2])
    return self._alignments[(c1,c2)]

  @property
  def chem_groups(self):
    """
    Intra-complex group of chemically identical (seq. id. > 95%) polypeptide
    chains as extracted from :attr:`ca_chains`. First chain in group is the one
    with the 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 = _GetChemGroups(self, 95)
      LogInfo('Chemically equivalent chain-groups in %s: %s' \
              % (self.GetName(), str(self._chem_groups)))
    return self._chem_groups

  def GetAngles(self, c1, c2):
    """Get Euler angles from superposition of chain *c1* with chain *c2*.
    Computed on first use based on :attr:`ca_chains` (cached).

    :param c1: Chain name for first chain to superpose.
    :type c1:  :class:`str`
    :param c2: Chain name for second chain to superpose.
    :type c2:  :class:`str`
    :return: 3 Euler angles (may contain nan if something fails).
    :rtype:  :class:`numpy.array`
    """
    if (c1,c2) not in self._angles:
      self._GetSuperposeData(c1, c2)
    return self._angles[(c1,c2)]

  def GetAxis(self, c1, c2):
    """Get axis of symmetry from superposition of chain *c1* with chain *c2*.
    Computed on first use based on :attr:`ca_chains` (cached).

    :param c1: Chain name for first chain to superpose.
    :type c1:  :class:`str`
    :param c2: Chain name for second chain to superpose.
    :type c2:  :class:`str`
    :return: Rotational axis (may contain nan if something fails).
    :rtype:  :class:`numpy.array`
    """
    if (c1,c2) not in self._axis:
      self._GetSuperposeData(c1, c2)
    return self._axis[(c1,c2)]

  @property
  def contacts(self):
    """
    Connectivity dictionary (**read/write**).
    As given by :func:`GetContacts` with *calpha_only* = False on :attr:`ent`.

    :getter: Computed on first use (cached)
    :setter: Uses :func:`FilterContacts` to ensure that we only keep contacts
             for chains in the cleaned entity.
    :type: See return type of :func:`GetContacts`
    """
    if self._contacts is None:
      self._contacts = GetContacts(self.ent, False)
    return self._contacts

  @contacts.setter
  def contacts(self, new_contacts):
    chain_names = set([ch.name for ch in self.ent.chains])
    self._contacts = FilterContacts(new_contacts, chain_names)
  
  @property
  def contacts_ca(self):
    """
    CA-only connectivity dictionary (**read/write**).
    Like :attr:`contacts` but with *calpha_only* = True in :func:`GetContacts`.
    """
    if self._contacts_ca is None:
      self._contacts_ca = GetContacts(self.ent, True)
    return self._contacts_ca
  
  @contacts_ca.setter
  def contacts_ca(self, new_contacts):
    chain_names = set([ch.name for ch in self.ent.chains])
    self._contacts_ca = FilterContacts(new_contacts, chain_names)

  ##############################################################################
  # Class internal helpers (anything that doesnt easily work without this class)
  ##############################################################################

  def _GetSuperposeData(self, c1, c2):
    """Fill _angles and _axis from superposition of CA chains of c1 and c2."""
    # get aligned views (must contain identical numbers of atoms!)
    aln = self.GetAlignment(c1, c2)
    if not aln:
      # fallback for non-aligned stuff (nan)
      self._angles[(c1,c2)] = np.empty(3) * np.nan
      self._axis[(c1,c2)] = np.empty(3) * np.nan
      return
    v1, v2 = seq.ViewsFromAlignment(aln)
    if v1.atom_count < 3:
      # fallback for non-aligned stuff (nan)
      self._angles[(c1,c2)] = np.empty(3) * np.nan
      self._axis[(c1,c2)] = np.empty(3) * np.nan
      return
    # get transformation
    sup_res = mol.alg.SuperposeSVD(v1, v2, apply_transform=False)
    Rt = sup_res.transformation
    # extract angles
    a,b,c = _GetAngles(Rt)
    self._angles[(c1,c2)] = np.asarray([a,b,c])
    # extract axis of symmetry
    R3 = geom.Rotation3(Rt.ExtractRotation())
    self._axis[(c1,c2)] = np.asarray(R3.GetRotationAxis().data)

###############################################################################
# Contacts computations
###############################################################################

def FilterContacts(contacts, chain_names):
  """Filter contacts to contain only contacts for chains in *chain_names*.

  :param contacts: Connectivity dictionary as produced by :func:`GetContacts`.
  :type contacts:  :class:`dict`
  :param chain_names: Chain names to keep.
  :type chain_names:  :class:`list` or (better) :class:`set`
  :return: New connectivity dictionary (format as in :func:`GetContacts`)
  :rtype:  :class:`dict`
  """
  # create new dict with subset
  filtered_contacts = dict()
  for c1 in contacts:
    if c1 in chain_names:
      new_contacts = dict()
      for c2 in contacts[c1]:
        if c2 in chain_names:
          new_contacts[c2] = contacts[c1][c2]
      # avoid adding empty dicts
      if new_contacts:
        filtered_contacts[c1] = new_contacts
  return filtered_contacts

def GetContacts(entity, calpha_only, dist_thr=12.0):
  """Get inter-chain contacts of a macromolecular entity.

  Contacts are pairs of residues within a given distance belonging to different
  chains. They are stored once per pair and include the CA/CB-CA/CB distance.

  :param entity: An entity to check connectivity for.
  :type entity:  :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
  :param calpha_only: If True, we only consider CA-CA distances. Else, we use CB
                      unless the residue is a GLY.
  :type calpha_only:  :class:`bool`
  :param dist_thr: Maximal CA/CB-CA/CB distance to be considered in contact.
  :type dist_thr:  :class:`float`
  :return: A connectivity dictionary. A pair of residues with chain names
           *ch_name1* & *ch_name2* (*ch_name1* < *ch_name2*), residue numbers
           *res_num1* & *res_num2* and distance *dist* (<= *dist_thr*) are
           stored as *result[ch_name1][ch_name2][res_num1][res_num2]* = *dist*.
  :rtype:  :class:`dict`
  """
  # get ent copy to search on
  if calpha_only:
    ev = entity.Select("aname=CA")
  else:
    ev = entity.Select("(rname=GLY and aname=CA) or aname=CB")
  ent = mol.CreateEntityFromView(ev, True)
  # search all vs all
  contacts = dict()
  for atom in ent.atoms:
    ch_name1 = atom.chain.name
    res_num1 = atom.residue.number.num
    close_atoms = ent.FindWithin(atom.pos, dist_thr)
    for close_atom in close_atoms:
      ch_name2 = close_atom.chain.name
      if ch_name2 > ch_name1:
        res_num2 = close_atom.residue.number.num
        dist = geom.Distance(atom.pos, close_atom.pos)
        # add to contacts
        if ch_name1 not in contacts:
          contacts[ch_name1] = dict()
        if ch_name2 not in contacts[ch_name1]:
          contacts[ch_name1][ch_name2] = dict()
        if res_num1 not in contacts[ch_name1][ch_name2]:
          contacts[ch_name1][ch_name2][res_num1] = dict()
        contacts[ch_name1][ch_name2][res_num1][res_num2] = round(dist, 3)
  # DONE
  return contacts


###############################################################################
# HELPERS
###############################################################################

# general

def _AlignAtomSeqs(seq_1, seq_2):
  """
  :type seq_1: :class:`ost.seq.SequenceHandle`
  :type seq_2: :class:`ost.seq.SequenceHandle`
  :return: Alignment of two sequences using a global aignment. Views attached
           to the input sequences will remain attached in the aln.
  :rtype:  :class:`~ost.seq.AlignmentHandle` or None if it failed.
  """
  # NOTE: If the two sequence have a greatly different length
  #       a local alignment could be a better choice...
  aln = None
  alns = seq.alg.GlobalAlign(seq_1, seq_2, seq.alg.BLOSUM62)
  if alns: aln = alns[0]
  if not aln:
    LogWarning('Failed to align %s to %s' % (seq_1.name, seq_2.name))
    LogWarning('%s:  %s' % (seq_1.name, seq_1.string))
    LogWarning('%s:  %s' % (seq_2.name, seq_2.string))
  return aln

# QS entity

def _CleanInputEntity(ent):
  """
  :param ent: The OST entity to be cleaned.
  :type ent:  :class:`EntityHandle` or :class:`EntityView`
  :return: A tuple of 3 items: :attr:`QSscoreEntity.ent`,
                               :attr:`QSscoreEntity.removed_chains`,
                               :attr:`QSscoreEntity.calpha_only`
  """
  # find chains to remove
  removed_chains = []
  for ch in ent.chains:
    # we remove chains if they are small-peptides or if the contain no aa
    # or if they contain only unknown or modified residues
    if    ch.name in ['-', '_'] \
       or ch.residue_count < 20 \
       or not any(r.chem_type.IsAminoAcid() for r in ch.residues) \
       or not (set(r.one_letter_code for r in ch.residues) - {'?', 'X'}):
      removed_chains.append(ch.name)

  # remove them from *ent*
  if removed_chains:
    chain_set = set()
    for ch_name in removed_chains:
      if ch_name in ['-', '_', ' ']:
        chain_set.add('"%c"' % ch_name)
      else:
        chain_set.add(ch_name)
    view = ent.Select('cname!=%s' % ','.join(chain_set))
    ent_new = mol.CreateEntityFromView(view, True)
    ent_new.SetName(ent.GetName())
  else:
    ent_new = ent

  # check if CA only
  calpha_only = False
  if ent_new.Select('aname=CB').atom_count == 0:
    LogInfo('Structure %s is a CA only structure!' % ent_new.GetName())
    calpha_only = True

  # report and return
  if removed_chains:
    LogInfo('Chains removed from %s: %s' \
            % (ent_new.GetName(), ''.join(removed_chains)))
  LogInfo('Chains in %s: %s' \
          % (ent_new.GetName(), ''.join([c.name for c in ent_new.chains])))
  return ent_new, removed_chains, calpha_only

def _GetCAOnlyEntity(ent):
  """
  :param ent: Entity to process
  :type ent:  :class:`EntityHandle` or :class:`EntityView`
  :return: New entity with only CA and only one atom per residue
           (see :attr:`QSscoreEntity.ca_entity`)
  """
  # cook up CA only view (diff from Select = guaranteed 1 atom per residue)
  ca_view = ent.CreateEmptyView()
  # add chain by chain
  for res in ent.residues:
    ca_atom = res.FindAtom("CA")
    if ca_atom.IsValid():
      ca_view.AddAtom(ca_atom)
  # finalize
  return mol.CreateEntityFromView(ca_view, False)

def _GetChemGroups(qs_ent, seqid_thr=95.):
  """
  :return: Intra-complex group of chemically identical polypeptide chains
           (see :attr:`QSscoreEntity.chem_groups`)

  :param qs_ent: Entity to process
  :type qs_ent:  :class:`QSscoreEntity`
  :param seqid_thr: Threshold used to decide when two chains are identical.
                    95 percent tolerates the few mutations crystallographers
                    like to do.
  :type seqid_thr:  :class:`float`
  """
  # get data from qs_ent
  ca_chains = qs_ent.ca_chains
  chain_names = sorted(ca_chains.keys())
  # get pairs of identical chains
  # NOTE: this scales quadratically with number of chains and may be optimized
  #       -> one could merge it with "merge transitive pairs" below...
  id_seqs = []
  for ch_1, ch_2 in itertools.combinations(chain_names, 2):
    aln = qs_ent.GetAlignment(ch_1, ch_2)
    if aln and seq.alg.SequenceIdentity(aln) > seqid_thr:
      id_seqs.append((ch_1, ch_2))
  # trivial case: no matching pairs
  if not id_seqs:
    return [[name] for name in chain_names]

  # merge transitive pairs
  groups = []
  for ch_1, ch_2 in id_seqs:
    found = False
    for g in groups:
      if ch_1 in g or ch_2 in g:
        found = True
        g.add(ch_1)
        g.add(ch_2)
    if not found:
      groups.append(set([ch_1, ch_2]))
  # sort internally based on sequence length
  chem_groups = []
  for g in groups:
    ranked_g = sorted([(-ca_chains[ch].length, ch) for ch in g])
    chem_groups.append([ch for _,ch in ranked_g])
  # add other dissimilar chains
  for ch in chain_names:
    if not any(ch in g for g in chem_groups):
      chem_groups.append([ch])
  
  return chem_groups

def _GetAngles(Rt):
  """Computes the Euler angles given a transformation matrix.

  :param Rt: Rt operator.
  :type Rt:  :class:`ost.geom.Mat4`
  :return: A :class:`tuple` of angles for each axis (x,y,z)
  """
  rot = np.asmatrix(Rt.ExtractRotation().data).reshape(3,3)
  tx = np.arctan2(rot[2,1], rot[2,2])
  if tx < 0:
    tx += 2*np.pi
  ty = np.arctan2(rot[2,0], np.sqrt(rot[2,1]**2 + rot[2,2]**2))
  if ty < 0:
    ty += 2*np.pi
  tz = np.arctan2(rot[1,0], rot[0,0])
  if tz < 0:
    tz += 2*np.pi
  return tx,ty,tz

# QS scorer

def _GetChemGroupsMapping(qs_ent_1, qs_ent_2):
  """
  :return: Inter-complex mapping of chemical groups
           (see :attr:`QSscorer.chem_mapping`)

  :param qs_ent_1: See :attr:`QSscorer.qs_ent_1`
  :param qs_ent_2: See :attr:`QSscorer.qs_ent_2`
  """
  # get chem. groups and unique representative
  chem_groups_1 = qs_ent_1.chem_groups
  chem_groups_2 = qs_ent_2.chem_groups
  repr_chains_1 = {x[0]: tuple(x) for x in chem_groups_1}
  repr_chains_2 = {x[0]: tuple(x) for x in chem_groups_2}

  # if entities do not have different number of unique chains, we get the
  # mapping for the smaller set
  swapped = False
  if len(repr_chains_2) < len(repr_chains_1):
    repr_chains_1, repr_chains_2 = repr_chains_2, repr_chains_1
    qs_ent_1, qs_ent_2 = qs_ent_2, qs_ent_1
    swapped = True

  # find the closest to each chain between the two entities
  # NOTE: this may still be sensible to orthology problem
  # -> currently we use a global alignment and seq. id. to rank pairs
  # -> we also tried local alignments and weighting the seq. id. by the
  #    coverages of the alignments (gapless string in aln. / seq. length)
  #    but global aln performed better...
  chain_pairs = []
  ca_chains_1 = qs_ent_1.ca_chains
  ca_chains_2 = qs_ent_2.ca_chains
  for ch_1 in repr_chains_1.keys():
    for ch_2 in repr_chains_2.keys():
      aln = _AlignAtomSeqs(ca_chains_1[ch_1], ca_chains_2[ch_2])
      if aln:
        chains_seqid = seq.alg.SequenceIdentity(aln)
        LogVerbose('Sequence identity', ch_1, ch_2, 'seqid=%.2f' % chains_seqid)
        chain_pairs.append((chains_seqid, ch_1, ch_2))

  # get top matching groups first
  chain_pairs = sorted(chain_pairs, reverse=True)
  chem_mapping = {}
  for _, c1, c2 in chain_pairs:
    skip = False
    for a,b in chem_mapping.iteritems():
      if repr_chains_1[c1] == a or repr_chains_2[c2] == b:
        skip = True
        break
    if not skip:
      chem_mapping[repr_chains_1[c1]] = repr_chains_2[c2]
  if swapped:
    chem_mapping = {y: x for x, y in chem_mapping.iteritems()}
    qs_ent_1, qs_ent_2 = qs_ent_2, qs_ent_1

  # notify chains without partner
  mapped_1 = set([i for s in chem_mapping.keys() for i in s])
  chains_1 = set([c.name for c in qs_ent_1.ent.chains])
  if chains_1 - mapped_1:
    LogWarning('Unmapped Chains in %s: %s'
               % (qs_ent_1.GetName(), ','.join(list(chains_1 - mapped_1))))

  mapped_2 = set([i for s in chem_mapping.values() for i in s])
  chains_2 = set([c.name for c in qs_ent_2.ent.chains])
  if chains_2 - mapped_2:
    LogWarning('Unmapped Chains in %s: %s'
               % (qs_ent_2.GetName(), ','.join(list(chains_2 - mapped_2))))
  
  # check if we have any chains left
  LogInfo('Chemical chain-groups mapping: ' + str(chem_mapping))
  if len(mapped_1) < 2 or len(mapped_2) < 2:
    raise QSscoreError('Less than 2 chains left in chem_mapping.')
  return chem_mapping

def _SelectFew(l, max_elements):
  """Return l or copy of l with at most *max_elements* entries."""
  n_el = len(l)
  if n_el <= max_elements:
    return l
  else:
    # cheap integer ceiling (-1 to ensure that x*max_elements gets d_el = x)
    d_el = ((n_el - 1) // max_elements) + 1
    new_l = list()
    for i in range(0, n_el, d_el):
      new_l.append(l[i])
    return new_l

def _GetAlignedResidues(qs_ent_1, qs_ent_2, chem_mapping, max_ca_per_chain,
                        clustalw_bin):
  """
  :return: Tuple of two :class:`~ost.mol.EntityView` objects containing subsets
           of *qs_ent_1* and *qs_ent_2*. Two entities are later created from
           those views (see :attr:`QSscorer.ent_to_cm_1` and
           :attr:`QSscorer.ent_to_cm_2`)

  :param qs_ent_1: See :attr:`QSscorer.qs_ent_1`
  :param qs_ent_2: See :attr:`QSscorer.qs_ent_2`
  :param chem_mapping: See :attr:`QSscorer.chem_mapping`
  :param max_ca_per_chain: See :attr:`QSscorer.max_ca_per_chain_for_cm`
  """
  # resulting views into CA entities using CA chain sequences
  ent_view_1 = qs_ent_1.ca_entity.CreateEmptyView()
  ent_view_2 = qs_ent_2.ca_entity.CreateEmptyView()
  ca_chains_1 = qs_ent_1.ca_chains
  ca_chains_2 = qs_ent_2.ca_chains
  # go through all mapped chemical groups
  for group_1, group_2 in chem_mapping.iteritems():
    # MSA with ClustalW
    seq_list = seq.CreateSequenceList()
    # keep sequence-name (must be unique) to view mapping
    seq_to_empty_view = dict()
    for ch in group_1:
      sequence = ca_chains_1[ch].Copy()
      sequence.name = qs_ent_1.GetName() + '.' + ch
      seq_to_empty_view[sequence.name] = ent_view_1
      seq_list.AddSequence(sequence)
    for ch in group_2:
      sequence = ca_chains_2[ch].Copy()
      sequence.name = qs_ent_2.GetName() + '.' + ch
      seq_to_empty_view[sequence.name] = ent_view_2
      seq_list.AddSequence(sequence)
    alnc = ClustalW(seq_list, clustalw=clustalw_bin)

    # collect aligned residues (list of lists of sequence_count valid residues)
    residue_lists = list()
    sequence_count = alnc.sequence_count
    for col in alnc:
      residues = [col.GetResidue(ir) for ir in range(sequence_count)]
      if all([r.IsValid() for r in residues]):
        residue_lists.append(residues)
    # check number of aligned residues
    if len(residue_lists) < 5:
      raise QSscoreError('Not enough aligned residues.')
    elif max_ca_per_chain > 0:
      residue_lists = _SelectFew(residue_lists, max_ca_per_chain)
    # check what view is for which residue
    res_views = [seq_to_empty_view[alnc.sequences[ir].name] \
                 for ir in range(sequence_count)]
    # add to view (note: only 1 CA atom per residue in here)
    for res_list in residue_lists:
      for res, view in zip(res_list, res_views):
        view.AddResidue(res, mol.ViewAddFlag.INCLUDE_ATOMS)
  # done
  return ent_view_1, ent_view_2


def _FindSymmetry(qs_ent_1, qs_ent_2, ent_to_cm_1, ent_to_cm_2, chem_mapping):
  """
  :return: A pair of comparable symmetry groups (for :attr:`QSscorer.symm_1`
           and :attr:`QSscorer.symm_2`) between the two structures.
           Empty lists if no symmetry identified.
  
  :param qs_ent_1: See :attr:`QSscorer.qs_ent_1`
  :param qs_ent_2: See :attr:`QSscorer.qs_ent_2`
  :param ent_to_cm_1: See :attr:`QSscorer.ent_to_cm_1`
  :param ent_to_cm_2: See :attr:`QSscorer.ent_to_cm_2`
  :param chem_mapping: See :attr:`QSscorer.chem_mapping`
  """
  LogInfo('Identifying Symmetry Groups...')

  # get possible symmetry groups
  symm_subg_1 = _GetSymmetrySubgroups(qs_ent_1, ent_to_cm_1,
                                      chem_mapping.keys())
  symm_subg_2 = _GetSymmetrySubgroups(qs_ent_2, ent_to_cm_2,
                                      chem_mapping.values())

  # check combination of groups
  LogInfo('Selecting Symmetry Groups...')
  # check how many mappings are to be checked for each pair of symmetry groups
  # -> lower number here will speed up worst-case runtime of _GetChainMapping
  # NOTE: this is tailored to speed up brute force all vs all mapping test
  #       for preferred _CheckClosedSymmetry this is suboptimal!
  best_symm = []
  for symm_1, symm_2 in itertools.product(symm_subg_1, symm_subg_2):
    s1 = symm_1[0]
    s2 = symm_2[0]
    if len(s1) != len(s2):
      LogDebug('Discarded', str(symm_1), str(symm_2),
               ': different size of symmetry groups')
      continue
    count = _CountSuperpositionsAndMappings(symm_1, symm_2, chem_mapping)
    nr_mapp = count['intra']['mappings'] + count['inter']['mappings']
    LogDebug('OK', str(symm_1), str(symm_2), ': %s mappings' % nr_mapp)
    best_symm.append((nr_mapp, symm_1, symm_2))

  for _, symm_1, symm_2 in sorted(best_symm):
    s1 = symm_1[0]
    s2 = symm_2[0]
    group_1 = ent_to_cm_1.Select('cname=%s' % ','.join(s1))
    group_2 = ent_to_cm_2.Select('cname=%s' % ','.join(s2))
    # check if by superposing a pair of chains within the symmetry group to
    # superpose all chains within the symmetry group
    # -> if successful, the symmetry groups are compatible
    closed_symm = _CheckClosedSymmetry(group_1, group_2, symm_1, symm_2,
                                       chem_mapping, 6, 0.8, find_best=False)

    if closed_symm:
      return symm_1, symm_2

  # nothing found
  LogInfo('No coherent symmetry identified between structures')
  return [], []


def _GetChainMapping(ent_1, ent_2, symm_1, symm_2, chem_mapping):
  """
  :return: Mapping from *ent_1* to *ent_2* (see :attr:`QSscorer.chain_mapping`)  

  :param ent_1: See :attr:`QSscorer.ent_to_cm_1`
  :param ent_2: See :attr:`QSscorer.ent_to_cm_2`
  :param symm_1: See :attr:`QSscorer.symm_1`
  :param symm_2: See :attr:`QSscorer.symm_2`
  :param chem_mapping: See :attr:`QSscorer.chem_mapping`
  """
  LogInfo('Symmetry-groups used in %s: %s' % (ent_1.GetName(), str(symm_1)))
  LogInfo('Symmetry-groups used in %s: %s' % (ent_2.GetName(), str(symm_2)))

  # quick check for closed symmetries
  thresholds = [(    'strict', 4, 0.8),
                (  'tolerant', 6, 0.4),
                ('permissive', 8, 0.2)]
  for scheme, sup_thr, sup_fract in thresholds:
    # check if by superposing a pair of chains within the symmetry group to
    # superpose all chains of the two oligomers
    # -> this also returns the resulting mapping of chains
    mapping = _CheckClosedSymmetry(ent_1, ent_2, symm_1, symm_2,
                                   chem_mapping, sup_thr, sup_fract)
    if mapping:
      LogInfo('Closed Symmetry with %s parameters' % scheme)
      if scheme == 'permissive':
        LogWarning('Permissive thresholds used for overlap based mapping ' + \
                   'detection: check mapping manually: %s' % mapping)
      return mapping
  
  # NOTE that what follows below is sub-optimal:
  # - if the two structures don't fit at all, we may map chains rather randomly
  # - we may need a lot of operations to find globally "best" mapping
  # -> COST = O(N^2) * O(N!)
  #    (first = all pairwise chain-superpose, latter = all chain mappings)
  # -> a greedy chain mapping choice algo (pick best RMSD for each one-by-one)
  #    could reduce this to O(N^2) * O(N^2) but it would need to be validated
  # - one could also try some other heuristic based on center of atoms of chains
  #   -> at this point we get a bad mapping anyways...

  # if we get here, we will need to check many combinations of mappings
  # -> check how many mappings must be checked and limit
  count = _CountSuperpositionsAndMappings(symm_1, symm_2, chem_mapping)
  LogInfo('Intra Symmetry-group mappings to check: %s' \
          % count['intra']['mappings'])
  LogInfo('Inter Symmetry-group mappings to check: %s' \
          % count['inter']['mappings'])
  nr_mapp = count['intra']['mappings'] + count['inter']['mappings']
  if nr_mapp > 100000: # 322560 is octamer vs octamer
    raise QSscoreError('Too many possible mappings: %s' % nr_mapp)

  # to speed up the computations we cache chain views and RMSDs
  cached_rmsd = _CachedRMSD(ent_1, ent_2)

  # get INTRA symmetry group chain mapping
  check = 0
  intra_mappings = [] # list of (RMSD, c1, c2, mapping) tuples (best superpose)
  # limit chem mapping to reference symmetry group
  ref_symm_1 = symm_1[0]
  ref_symm_2 = symm_2[0]
  asu_chem_mapping = _LimitChemMapping(chem_mapping, ref_symm_1, ref_symm_2)
  # for each chemically identical group
  for g1, g2 in asu_chem_mapping.iteritems():
    # try to superpose all possible pairs
    for c1, c2 in itertools.product(g1, g2):
      # get superposition transformation
      LogVerbose('Superposing chains: %s, %s' % (c1, c2))
      res = cached_rmsd.GetSuperposition(c1, c2)
      # compute RMSD of possible mappings
      cur_mappings = [] # list of (RMSD, mapping) tuples
      for mapping in _ListPossibleMappings(c1, c2, asu_chem_mapping):
        check += 1
        chains_rmsd = cached_rmsd.GetMappedRMSD(mapping, res.transformation)
        cur_mappings.append((chains_rmsd, mapping))
        LogVerbose(str(mapping), str(chains_rmsd))
      best_rmsd, best_mapp = min(cur_mappings)
      intra_mappings.append((best_rmsd, c1, c2, best_mapp))
  # best chain-chain superposition defines the intra asu mapping
  _, intra_asu_c1, intra_asu_c2, intra_asu_mapp = min(intra_mappings)

  # if only one asu is present in any of the complexes, we're done
  if len(symm_1) == 1 or len(symm_2) == 1:
    mapping = intra_asu_mapp
  else:
    # the mapping is the element position within the asu chain list
    # -> this speed up a lot, assuming that the order of chain in asu groups
    #    is following the order of symmetry operations
    index_asu_c1 = ref_symm_1.index(intra_asu_c1)
    index_asu_c2 = ref_symm_2.index(intra_asu_c2)
    index_mapp = {ref_symm_1.index(s1): ref_symm_2.index(s2) \
                  for s1, s2 in intra_asu_mapp.iteritems()}
    LogInfo('Intra symmetry-group mapping: %s' % str(intra_asu_mapp))

    # get INTER symmetry group chain mapping
    # we take the first symmetry group from the complex having the most
    if len(symm_1) < len(symm_2):
      symm_combinations = list(itertools.product(symm_1, [symm_2[0]]))
    else:
      symm_combinations = list(itertools.product([symm_1[0]], symm_2))

    full_asu_mapp = {tuple(symm_1): tuple(symm_2)}
    full_mappings = [] # list of (RMSD, mapping) tuples
    for s1, s2 in symm_combinations:
      # get superposition transformation (we take best chain-pair in asu)
      LogVerbose('Superposing symmetry-group: %s, %s in %s, %s' \
                 % (s1[index_asu_c1], s2[index_asu_c2], s1, s2))
      res = cached_rmsd.GetSuperposition(s1[index_asu_c1], s2[index_asu_c2])
      # compute RMSD of possible mappings
      for asu_mapp in _ListPossibleMappings(s1, s2, full_asu_mapp):
        check += 1
        # need to extract full chain mapping (use indexing)
        mapping = {}
        for a, b in asu_mapp.iteritems():
          for id_a, id_b in index_mapp.iteritems():
            mapping[a[id_a]] = b[id_b]
        chains_rmsd = cached_rmsd.GetMappedRMSD(mapping, res.transformation)
        full_mappings.append((chains_rmsd, mapping))
        LogVerbose(str(mapping), str(chains_rmsd))

    if check != nr_mapp:
      LogWarning('Mapping number estimation was wrong') # sanity check

    # return best (lowest RMSD) mapping
    mapping = min(full_mappings)[1]

  LogWarning('Extensive search used for mapping detection (fallback). This ' + \
             'approach has known limitations. Check mapping manually: %s' \
             % mapping)
  return mapping


def _GetSymmetrySubgroups(qs_ent, ent, chem_groups):
  """
  Identify the symmetry (either cyclic C or dihedral D) of the protein and find
  all possible symmetry subgroups. This is done testing all combination of chain
  superposition and clustering them by the angles (D) or the axis (C) of the Rt
  operator.

  Groups of superposition which can fully reconstruct the structure are possible
  symmetry subgroups.

  :param qs_ent: Entity with cached angles and axis.
  :type qs_ent:  :class:`QSscoreEntity`
  :param ent: Entity from which to extract chains (perfect alignment assumed
              for superposition as in :attr:`QSscorer.ent_to_cm_1`).
  :type ent:  :class:`EntityHandle` or :class:`EntityView`
  :param chem_groups: List of tuples/lists of chain names in *ent*. Each list
                      contains all chains belonging to a chem. group (could be
                      from :attr:`QSscoreEntity.chem_groups` or from "keys()"
                      or "values()" of :attr:`QSscorer.chem_mapping`)

  :return: A list of possible symmetry subgroups (each in same format as
           :attr:`QSscorer.symm_1`). If no symmetry is found, we return a list
           with a single symmetry subgroup with a single group with all chains.
  """
  # get angles / axis for pairwise transformations within same chem. group
  angles = {}
  axis = {}
  for cnames in chem_groups:
    for c1, c2 in itertools.combinations(cnames, 2):
      cur_angles = qs_ent.GetAngles(c1, c2)
      if not np.any(np.isnan(cur_angles)):
        angles[(c1,c2)] = cur_angles
      cur_axis = qs_ent.GetAxis(c1, c2)
      if not np.any(np.isnan(cur_axis)):
        axis[(c1,c2)] = cur_axis

  # cluster superpositions angles at different thresholds
  symm_groups = []
  LogVerbose('Possible symmetry-groups for: %s' % ent.GetName())
  for angle_thr in np.arange(0.1, 1, 0.1):
    d_symm_groups = _GetDihedralSubgroups(ent, chem_groups, angles, angle_thr)
    if d_symm_groups:
      for group in d_symm_groups:
        if group not in symm_groups:
          symm_groups.append(group)
          LogVerbose('Dihedral: %s' % group)
      LogInfo('Symmetry threshold %.1f used for angles of %s' \
              % (angle_thr, ent.GetName()))
      break
  
  # cluster superpositions axis at different thresholds
  for axis_thr in np.arange(0.1, 1, 0.1):
    c_symm_groups = _GetCyclicSubgroups(ent, chem_groups, axis, axis_thr)
    if c_symm_groups:
      for group in c_symm_groups:
        if group not in symm_groups:
          symm_groups.append(group)
          LogVerbose('Cyclic: %s' % group)
      LogInfo('Symmetry threshold %.1f used for axis of %s' \
              % (axis_thr, ent.GetName()))
      break

  # fall back to single "group" containing all chains if none found
  if not symm_groups:
    LogInfo('No symmetry-group detected in %s' % ent.GetName())
    symm_groups = [[tuple([c for g in chem_groups for c in g])]]

  return symm_groups

def _GetDihedralSubgroups(ent, chem_groups, angles, angle_thr):
  """
  :return: Dihedral subgroups for :func:`_GetSymmetrySubgroups`
           (same return type as there). Empty list if fail.

  :param ent: See :func:`_GetSymmetrySubgroups`
  :param chem_groups: See :func:`_GetSymmetrySubgroups`
  :param angles: :class:`dict` (key = chain-pair-tuple, value = Euler angles)
  :param angle_thr: Euler angles distance threshold for clustering (float).
  """
  # cluster superpositions angles
  if len(angles) == 0:
    # nothing to be done here
    return []
  else:
    same_angles = _ClusterData(angles, angle_thr, _AngleArrayDistance)

  # get those which are non redundant and covering all chains
  symm_groups = []
  for clst in same_angles.values():
    group = clst.keys()
    if _ValidChainGroup(group, ent):
      if len(chem_groups) > 1:
        # if hetero, we want to group toghether different chains only
        symm_groups.append(zip(*group))
      else:
        # if homo, we also want pairs
        symm_groups.append(group)
        symm_groups.append(zip(*group))
  return symm_groups

def _GetCyclicSubgroups(ent, chem_groups, axis, axis_thr):
  """
  :return: Cyclic subgroups for :func:`_GetSymmetrySubgroups`
           (same return type as there). Empty list if fail.

  :param ent: See :func:`_GetSymmetrySubgroups`
  :param chem_groups: See :func:`_GetSymmetrySubgroups`
  :param angles: :class:`dict` (key = chain-pair-tuple, value = rotation axis)
  :param angle_thr: Axis distance threshold for clustering (float).
  """
  # cluster superpositions axis
  if len(axis) == 0:
    # nothing to be done here
    return []
  else:
    same_axis = _ClusterData(axis, axis_thr, _AxisDistance)

  # use to get grouping
  symm_groups = []
  for clst in same_axis.values():
    all_chain = [item for sublist in clst.keys() for item in sublist]
    if len(set(all_chain)) == ent.chain_count:
      # if we have an hetero we can exploit cyclic symmetry for grouping
      if len(chem_groups) > 1:
        ref_chem_group = chem_groups[0]
        # try two criteria for grouping
        cyclic_group_closest = []
        cyclic_group_iface = []
        for c1 in ref_chem_group:
          group_closest = [c1]
          group_iface = [c1]
          for chains in chem_groups[1:]:
            # center of atoms distance
            closest = _GetClosestChain(ent, c1, chains)
            # chain with bigger interface
            closest_iface = _GetClosestChainInterface(ent, c1, chains)
            group_closest.append(closest)
            group_iface.append(closest_iface)
          cyclic_group_closest.append(tuple(group_closest))
          cyclic_group_iface.append(tuple(group_iface))
        if _ValidChainGroup(cyclic_group_closest, ent):
          symm_groups.append(cyclic_group_closest)
        if _ValidChainGroup(cyclic_group_iface, ent):
          symm_groups.append(cyclic_group_iface)
      # if it is a homo, there's not much we can group
      else:
        symm_groups.append(chem_groups)
  return symm_groups

def _ClusterData(data, thr, metric):
  """Wrapper for fclusterdata to get dict of clusters.
  
  :param data: :class:`dict` (keys for ID, values used for clustering)
  :return: :class:`dict` {cluster_idx: {data-key: data-value}}
  """
  # special case length 1
  if len(data) == 1:
    return {0: {data.keys()[0]: data.values()[0]} }
  # do the clustering
  cluster_indices = fclusterdata(np.asarray(data.values()), thr,
                                 method='complete', criterion='distance',
                                 metric=metric)
  # fclusterdata output is cluster ids -> put into dict of clusters
  res = {}
  for cluster_idx, data_key in zip(cluster_indices, data.keys()):
    if not cluster_idx in res:
      res[cluster_idx] = {}
    res[cluster_idx][data_key] = data[data_key]
  return res

def _AngleArrayDistance(u, v):
  """
  :return: Average angular distance of two arrays of angles.
  :param u: Euler angles.
  :param v: Euler angles.
  """
  dist = []
  for a,b in zip(u,v):
    delta = abs(a-b)
    if delta > np.pi:
      delta = abs(2*np.pi - delta)
    dist.append(delta)
  return np.mean(dist)

def _AxisDistance(u, v):
  """
  :return: Euclidean distance between two rotation axes. Axes can point in
           either direction so we ensure to use the closer one.
  :param u: Rotation axis.
  :param v: Rotation axis.
  """
  # get axes as arrays (for convenience)
  u = np.asarray(u)
  v = np.asarray(v)
  # get shorter of two possible distances (using v or -v)
  dv1 = u - v
  dv2 = u + v
  d1 = np.dot(dv1, dv1)
  d2 = np.dot(dv2, dv2)
  return np.sqrt(min(d1, d2))

def _GetClosestChain(ent, ref_chain, chains):
  """
  :return: Chain closest to *ref_chain* based on center of atoms distance.
  :rtype:  :class:`str`
  :param ent: See :func:`_GetSymmetrySubgroups`
  :param ref_chain: We look for chains closest to this one
  :type ref_chain:  :class:`str`
  :param chains: We only consider these chains
  :type chains:  :class:`list` of :class:`str`
  """
  contacts = []
  ref_center = ent.FindChain(ref_chain).center_of_atoms
  for ch in chains:
    center = ent.FindChain(ch).center_of_atoms
    contacts.append((geom.Distance(ref_center, center), ch))
  closest_chain = min(contacts)[1]
  return closest_chain

def _GetClosestChainInterface(ent, ref_chain, chains):
  """
  :return: Chain with biggest interface (within 10 A) to *ref_chain*.
  :rtype:  :class:`str`
  :param ent: See :func:`_GetSymmetrySubgroups`
  :param ref_chain: We look for chains closest to this one
  :type ref_chain:  :class:`str`
  :param chains: We only consider these chains
  :type chains:  :class:`list` of :class:`str`
  """
  # NOTE: this is computed on a subset of the full biounit and might be
  # inaccurate. Also it could be extracted from QSscoreEntity.contacts.
  closest = []
  for ch in chains:
    iface_view = ent.Select('cname=%s and 10 <> [cname=%s]' % (ref_chain, ch))
    nr_res = iface_view.residue_count
    closest.append((nr_res, ch))
  closest_chain = max(closest)[1]
  return closest_chain

def _ValidChainGroup(group, ent):
  """
  :return: True, if *group* has unique chain names and as many chains as *ent*
  :rtype:  :class:`bool`
  :param group: Symmetry groups to check
  :type group:  Same as :attr:`QSscorer.symm_1`
  :param ent: See :func:`_GetSymmetrySubgroups`
  """
  all_chain = [item for sublist in group for item in sublist]
  if len(all_chain) == len(set(all_chain)) and\
     len(all_chain) == ent.chain_count:
    return True
  return False

def _LimitChemMapping(chem_mapping, limit_1, limit_2):
  """
  :return: Chem. mapping containing only chains in *limit_1* and *limit_2*
  :rtype:  Same as :attr:`QSscorer.chem_mapping`
  :param chem_mapping: See :attr:`QSscorer.chem_mapping`
  :param limit_1: Limits chain names in chem_mapping.keys()
  :type limit_1:  List/tuple of strings
  :param limit_2: Limits chain names in chem_mapping.values()
  :type limit_2:  List/tuple of strings
  """
  # use set intersection to keep only mapping for chains in limit_X
  limit_1_set = set(limit_1)
  limit_2_set = set(limit_2)
  asu_chem_mapping = dict()
  for key, value in chem_mapping.iteritems():
    new_key = tuple(set(key) & limit_1_set)
    if new_key:
      asu_chem_mapping[new_key] = tuple(set(value) & limit_2_set)
  return asu_chem_mapping


def _CountSuperpositionsAndMappings(symm_1, symm_2, chem_mapping):
  """
  :return: Dictionary of number of mappings and superpositions to be performed.
           Returned as *result[X][Y] = number* with X = "intra" or "inter" and
           Y = "mappings" or "superpositions". The idea is that for each
           pairwise superposition we check all possible mappings.
           We can check the combinations within (intra) a symmetry group and
           once established, we check the combinations between different (inter)
           symmetry groups.
  :param symm_1: See :attr:`QSscorer.symm_1`
  :param symm_2: See :attr:`QSscorer.symm_2`
  :param chem_mapping: See :attr:`QSscorer.chem_mapping`
  """
  # setup
  c = {}
  c['intra'] = {}
  c['inter'] = {}
  c['intra']['mappings'] = 0
  c['intra']['superpositions'] = 0
  c['inter']['mappings'] = 0
  c['inter']['superpositions'] = 0
  # intra ASU mappings within reference symmetry groups
  ref_symm_1 = symm_1[0]
  ref_symm_2 = symm_2[0]
  asu_chem_mapping = _LimitChemMapping(chem_mapping, ref_symm_1, ref_symm_2)
  for g1, g2 in asu_chem_mapping.iteritems():
    chain_superpositions = len(g1) * len(g2)
    c['intra']['superpositions'] += chain_superpositions
    map_per_sup = _PermutationOrCombinations(g1[0], g2[0], asu_chem_mapping)
    c['intra']['mappings'] += chain_superpositions * map_per_sup
  if len(symm_1) == 1 or len(symm_2) == 1:
    return c
  # inter ASU mappings
  asu_superposition = min(len(symm_1), len(symm_2))
  c['inter']['superpositions'] = asu_superposition
  asu = {tuple(symm_1): tuple(symm_2)}
  map_per_sup = _PermutationOrCombinations(ref_symm_1, ref_symm_2, asu)
  c['inter']['mappings'] = asu_superposition * map_per_sup
  return c

def _PermutationOrCombinations(sup1, sup2, chem_mapping):
  """Should match len(_ListPossibleMappings(sup1, sup2, chem_mapping))."""
  # remove superposed elements, put smallest complex as key
  chem_map = {}
  for a,b in chem_mapping.iteritems():
    new_a = tuple([x for x in a if x != sup1])
    new_b = tuple([x for x in b if x != sup2])
    chem_map[new_a] = new_b
  mapp_nr = 1.0
  for a, b in chem_map.iteritems():
    if len(a) == len(b):
      mapp_nr *= factorial(len(a))
    elif len(a) < len(b):
      mapp_nr *= binom(len(b), len(a))
    elif len(a) > len(b):
      mapp_nr *= binom(len(a), len(b))
  return int(mapp_nr)

def _ListPossibleMappings(sup1, sup2, chem_mapping):
  """
  Return a flat list of all possible mappings given *chem_mapping* and keeping
  mapping of *sup1* and *sup2* fixed. For instance if elements are chain names
  this is all the mappings to check for a given superposition.

  Elements in first complex are defined by *chem_mapping.keys()* (list of list
  of elements) and elements in second complex by *chem_mapping.values()*. If
  complexes don't have same number of elements, we map only elements for the one
  with less. Also mapping is only between elements of mapped groups according to
  *chem_mapping*.
           
  :rtype:  :class:`list` of :class:`dict` (key = element in chem_mapping-key,
           value = element in chem_mapping-value)
  :param sup1: Element for which mapping is fixed.
  :type sup1:  Like element in chem_mapping-key
  :param sup2: Element for which mapping is fixed.
  :type sup2:  Like element in chem_mapping-value
  :param chem_mapping: Defines mapping between groups of elements (e.g. result
                       of :func:`_LimitChemMapping`).
  :type chem_mapping:  :class:`dict` with key / value = :class:`tuple`

  :raises: :class:`QSscoreError` if reference complex (first one or one with
           less elements) has more elements for any given mapped group.
  """
  # find smallest complex
  chain1 = [i for s in chem_mapping.keys() for i in s]
  chain2 = [i for s in chem_mapping.values() for i in s]
  swap = False
  if len(chain1) > len(chain2):
    swap = True
  # remove superposed elements, put smallest complex as key
  chem_map = {}
  for a, b in chem_mapping.iteritems():
    new_a = tuple([x for x in a if x != sup1])
    new_b = tuple([x for x in b if x != sup2])
    if swap:
      chem_map[new_b] = new_a
    else:
      chem_map[new_a] = new_b
  # generate permutations or combinations of chemically
  # equivalent chains
  chem_perm = []
  chem_ref = []
  for a, b in chem_map.iteritems():
    if len(a) == len(b):
      chem_perm.append(list(itertools.permutations(b)))
      chem_ref.append(a)
    elif len(a) < len(b):
      chem_perm.append(list(itertools.combinations(b, len(a))))
      chem_ref.append(a)
    else:
      raise QSscoreError('Impossible to define reference group: %s' \
                         % str(chem_map))
  # save the list of possible mappings
  mappings = []
  flat_ref = [i for s in chem_ref for i in s]
  for perm in itertools.product(*chem_perm):
    flat_perm = [i for s in perm for i in s]
    d = {c1: c2 for c1, c2 in zip(flat_ref, flat_perm)}
    if swap:
      d = {v: k for k, v in d.items()}
    d.update({sup1: sup2})
    mappings.append(d)
  return mappings


def _CheckClosedSymmetry(ent_1, ent_2, symm_1, symm_2, chem_mapping,
                         sup_thr=4, sup_fract=0.8, find_best=True):
  """
  Quick check if we can superpose two chains and get a mapping for all other
  chains using the same transformation. The mapping is defined by sufficient
  overlap of the transformed chain of *ent_1* onto another chain in *ent_2*.

  :param ent_1: Entity to map to *ent_2* (perfect alignment assumed between
                chains of same chem. group as in :attr:`QSscorer.ent_to_cm_1`).
                Views are ok but only to select full chains!
  :param ent_2: Entity to map to (perfect alignment assumed between
                chains of same chem. group as in :attr:`QSscorer.ent_to_cm_2`).
                Views are ok but only to select full chains!
  :param symm_1: Symmetry groups to use. We only superpose chains within
                 reference symmetry group of *symm_1* and *symm_2*.
                 See :attr:`QSscorer.symm_1`
  :param symm_2: See :attr:`QSscorer.symm_2`
  :param chem_mapping: See :attr:`QSscorer.chem_mapping`.
                       All chains in *ent_1* / *ent_2* must be listed here!
  :param sup_thr: Distance around transformed chain in *ent_1* to check for
                  overlap.
  :type sup_thr:  :class:`float`
  :param sup_fract: Fraction of atoms in chain of *ent_2* that must be
                    overlapped for overlap to be sufficient.
  :type sup_fract:  :class:`float`
  :param find_best: If True, we look for best mapping according to
                    :func:`_ChainRMSD`. Otherwise, we return first suitable
                    mapping.
  :type find_best:  :class:`bool`

  :return: Mapping from *ent_1* to *ent_2* or None if none found. Mapping, if
           found, is as in :attr:`QSscorer.chain_mapping`.
  :rtype:  :class:`dict` (:class:`str` / :class:`str`)
  """
  # limit chem mapping to reference symmetry group
  ref_symm_1 = symm_1[0]
  ref_symm_2 = symm_2[0]
  asu_chem_mapping = _LimitChemMapping(chem_mapping, ref_symm_1, ref_symm_2)
  # for each chemically identical group
  rmsd_mappings = []
  for g1, g2 in asu_chem_mapping.iteritems():
    # try to superpose all possible pairs
    # -> note that some chain-chain combinations may work better than others
    #    to superpose the full oligomer (e.g. if some chains are open/closed)
    for c1, c2 in itertools.product(g1, g2):
      # get superposition transformation
      chain_1 = ent_1.Select('cname=%s' % c1)
      chain_2 = ent_2.Select('cname=%s' % c2)
      res = mol.alg.SuperposeSVD(chain_1, chain_2, apply_transform=False)
      # look for overlaps
      mapping = _GetSuperpositionMapping(ent_1, ent_2, chem_mapping,
                                         res.transformation, sup_thr,
                                         sup_fract)
      if not mapping:
        continue
      # early abort if we only want the first one
      if not find_best:
        return mapping
      else:
        # get RMSD for sorting
        rmsd = _GetMappedRMSD(ent_1, ent_2, mapping, res.transformation)
        rmsd_mappings.append((rmsd, mapping))
  # return best mapping
  if rmsd_mappings:
    return min(rmsd_mappings)[1]
  else:
    return None

def _GetSuperpositionMapping(ent_1, ent_2, chem_mapping, transformation,
                             sup_thr, sup_fract):
  """
  :return: Dict with chain mapping from *ent_1* to *ent_2* or None if failed
           (see :func:`_CheckClosedSymmetry`).
  :param ent_1: See :func:`_CheckClosedSymmetry`
  :param ent_2: See :func:`_CheckClosedSymmetry`
  :param chem_mapping: See :func:`_CheckClosedSymmetry`
  :param transformation: Superposition transformation to be applied to *ent_1*.
  :param sup_thr: See :func:`_CheckClosedSymmetry`
  :param sup_fract: See :func:`_CheckClosedSymmetry`
  """
  # NOTE: this seems to be the comput. most expensive part in QS scoring
  # -> it could be moved to C++ for speed-up
  # -> the next expensive bits are ClustalW and GetContacts btw...

  # swap if needed (ent_1 must be shorter or same)
  if ent_1.chain_count > ent_2.chain_count:
    swap = True
    ent_1, ent_2 = ent_2, ent_1
    transformation = geom.Invert(transformation)
    chem_pairs = zip(chem_mapping.values(), chem_mapping.keys())
  else:
    swap = False
    chem_pairs = chem_mapping.iteritems()
  # sanity check
  if ent_1.chain_count == 0:
    return None
  # extract valid partners for each chain in ent_1
  chem_partners = dict()
  for cg1, cg2 in chem_pairs:
    for ch in cg1:
      chem_partners[ch] = set(cg2)
  # get mapping for each chain in ent_1
  mapping = dict()
  mapped_chains = set()
  for ch_1 in ent_1.chains:
    # get all neighbors split by chain (NOTE: this muight be moved to C++)
    ch_atoms = {ch_2.name: set() for ch_2 in ent_2.chains}
    for a_1 in ch_1.handle.atoms:
      transformed_pos = geom.Vec3(transformation * geom.Vec4(a_1.pos))
      a_within = ent_2.FindWithin(transformed_pos, sup_thr)
      for a_2 in a_within:
        ch_atoms[a_2.chain.name].add(a_2.hash_code)
    # get one with most atoms in overlap
    max_num, max_name = max((len(atoms), name)
                            for name, atoms in ch_atoms.iteritems())
    # early abort if overlap fraction not good enough
    ch_2 = ent_2.FindChain(max_name)
    if max_num == 0 or max_num / float(ch_2.handle.atom_count) < sup_fract:
      return None
    # early abort if mapped to chain of different chem. group
    ch_1_name = ch_1.name
    if ch_1_name not in chem_partners:
      raise RuntimeError("Chem. mapping doesn't contain all chains!")
    if max_name not in chem_partners[ch_1_name]:
      return None
    # early abort if multiple ones mapped to same chain
    if max_name in mapped_chains:
      return None
    # set mapping
    mapping[ch_1_name] = max_name
    mapped_chains.add(max_name)
  # unswap if needed and return
  if swap:
    mapping = {v: k for k, v in mapping.iteritems()}
  return mapping

def _GetMappedRMSD(ent_1, ent_2, chain_mapping, transformation):
  """
  :return: RMSD between complexes considering chain mapping.
  :param ent_1: Entity mapped to *ent_2* (perfect alignment assumed between
                mapped chains as in :attr:`QSscorer.ent_to_cm_1`).
  :param ent_2: Entity which was mapped to (perfect alignment assumed between
                mapped chains as in :attr:`QSscorer.ent_to_cm_2`).
  :param chain_mapping: See :attr:`QSscorer.chain_mapping`
  :param transformation: Superposition transformation to be applied to *ent_1*.
  """
  # collect RMSDs and atom counts chain by chain and combine afterwards
  rmsds = []
  atoms = []
  for c1, c2 in chain_mapping.iteritems():
    # get views and atom counts
    chain_1 = ent_1.Select('cname=%s' % c1)
    chain_2 = ent_2.Select('cname=%s' % c2)
    atom_count = chain_1.atom_count
    if atom_count != chain_2.atom_count:
      raise RuntimeError('Chains in _GetMappedRMSD must be perfectly aligned!')
    # get RMSD
    rmsd = mol.alg.CalculateRMSD(chain_1, chain_2, transformation)
    # update lists
    rmsds.append(rmsd)
    atoms.append(float(atom_count))
  # combine (reminder: RMSD = sqrt(sum(atom_distance^2)/atom_count))
  rmsd = np.sqrt( sum([a * r**2 for a, r in zip(atoms, rmsds)]) / sum(atoms) )
  return rmsd

class _CachedRMSD:
  """Helper object for repetitive RMSD calculations.
  Meant to speed up :func:`_GetChainMapping` but could also be used to replace
  :func:`_GetMappedRMSD` in :func:`_CheckClosedSymmetry`.

  :param ent_1: See :attr:`QSscorer.ent_to_cm_1`
  :param ent_2: See :attr:`QSscorer.ent_to_cm_2`
  """
  def __init__(self, ent_1, ent_2):
    # initialize caches and keep entities
    self.ent_1 = ent_1
    self.ent_2 = ent_2
    self._chain_views_1 = dict()
    self._chain_views_2 = dict()
    self._pair_rmsd = dict()  # key = (c1,c2), value = (atom_count,rmsd)

  def GetChainView1(self, cname):
    """Get cached view on chain *cname* for :attr:`ent_1`."""
    if cname not in self._chain_views_1:
      self._chain_views_1[cname] = self.ent_1.Select('cname=%s' % cname)
    return self._chain_views_1[cname]

  def GetChainView2(self, cname):
    """Get cached view on chain *cname* for :attr:`ent_2`."""
    if cname not in self._chain_views_2:
      self._chain_views_2[cname] = self.ent_2.Select('cname=%s' % cname)
    return self._chain_views_2[cname]

  def GetSuperposition(self, c1, c2):
    """Get superposition result (no change in entities!) for *c1* to *c2*.
    This invalidates cached RMSD results used in :func:`GetMappedRMSD`.

    :param c1: chain name for :attr:`ent_1`.
    :param c2: chain name for :attr:`ent_2`.
    """
    # invalidate _pair_rmsd
    self._pair_rmsd = dict()
    # superpose
    chain_1 = self.GetChainView1(c1)
    chain_2 = self.GetChainView2(c2)
    if chain_1.atom_count != chain_2.atom_count:
      raise RuntimeError('Chains in GetSuperposition not perfectly aligned!')
    return mol.alg.SuperposeSVD(chain_1, chain_2, apply_transform=False)

  def GetMappedRMSD(self, chain_mapping, transformation):
    """
    :return: RMSD between complexes considering chain mapping.
    :param chain_mapping: See :attr:`QSscorer.chain_mapping`.
    :param transformation: Superposition transformation (e.g. res.transformation
                           for res = :func:`GetSuperposition`).
    """
    # collect RMSDs and atom counts chain by chain and combine afterwards
    rmsds = []
    atoms = []
    for c1, c2 in chain_mapping.iteritems():
      # cached?
      if (c1, c2) in self._pair_rmsd:
        atom_count, rmsd = self._pair_rmsd[(c1, c2)]
      else:
        # get views and atom counts
        chain_1 = self.GetChainView1(c1)
        chain_2 = self.GetChainView2(c2)
        atom_count = chain_1.atom_count
        if atom_count != chain_2.atom_count:
          raise RuntimeError('Chains in GetMappedRMSD not perfectly aligned!')
        # get RMSD
        rmsd = mol.alg.CalculateRMSD(chain_1, chain_2, transformation)
        self._pair_rmsd[(c1, c2)] = (atom_count, rmsd)
      # update lists
      rmsds.append(rmsd)
      atoms.append(float(atom_count))
    # combine (reminder: RMSD = sqrt(sum(atom_distance^2)/atom_count))
    rmsd = np.sqrt( sum([a * r**2 for a, r in zip(atoms, rmsds)]) / sum(atoms) )
    return rmsd


def _CleanUserSymmetry(symm, ent):
  """Clean-up user provided symmetry.

  :param symm: Symmetry group as in :attr:`QSscorer.symm_1`
  :param ent: Entity from which to extract chain names
  :type ent:  :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
  :return: New symmetry group limited to chains in sub-structure *ent*
           (see :attr:`QSscorer.symm_1`). Empty list if invalid symmetry.
  """
  # restrict symm to only contain chains in ent
  chain_names = set([ch.name for ch in ent.chains])
  new_symm = list()
  for symm_group in symm:
    new_group = tuple(ch for ch in symm_group if ch in chain_names)
    if new_group:
      new_symm.append(new_group)
  # report it
  if new_symm != symm:
    LogInfo("Cleaned user symmetry for %s: %s" % (ent.GetName(), new_symm))
  # do all groups have same length?
  lengths = [len(symm_group) for symm_group in new_symm]
  if lengths[1:] != lengths[:-1]:
    LogWarning('User symmetry group of different sizes for %s. Ignoring it!' \
               % (ent.GetName()))
    return []
  # do we cover all chains?
  s_chain_names = set([ch for symm_group in new_symm for ch in symm_group])
  if len(s_chain_names) != sum(lengths):
    LogWarning('User symmetry group for %s has duplicate chains. Ignoring it!' \
               % (ent.GetName()))
    return []
  if s_chain_names != chain_names:
    LogWarning('User symmetry group for %s is missing chains. Ignoring it!' \
               % (ent.GetName()))
    return []
  # ok all good
  return new_symm

def _AreValidSymmetries(symm_1, symm_2):
  """Check symmetry pair for major problems.

  :return: False if any of the two is empty or if they're compatible in size.
  :rtype:  :class:`bool`
  """
  if not symm_1 or not symm_2:
    return False
  if len(symm_1) != 1 or len(symm_2) != 1:
    if not len(symm_1[0]) == len(symm_2[0]):
      LogWarning('Symmetry groups of different sizes. Ignoring them!')
      return False
  return True

def _GetMappedAlignments(ent_1, ent_2, chain_mapping):
  """
  :return: Alignments of 2 structures given chain mapping
           (see :attr:`QSscorer.alignments`).
  :param ent_1: Entity containing all chains in *chain_mapping.keys()*.
                Views to this entity attached to first sequence of each aln.
  :param ent_2: Entity containing all chains in *chain_mapping.values()*.
                Views to this entity attached to second sequence of each aln.
  :param chain_mapping: See :attr:`QSscorer.chain_mapping`
  """
  alns = []
  for ch_1_name in sorted(chain_mapping):
    # get both sequences incl. attached view
    ch_1 = ent_1.FindChain(ch_1_name)
    seq_1 = seq.SequenceFromChain(ch_1.name, ch_1)
    ch_2 = ent_2.FindChain(chain_mapping[ch_1_name])
    seq_2 = seq.SequenceFromChain(ch_2.name, ch_2)
    # align them
    aln = _AlignAtomSeqs(seq_1, seq_2)
    if aln: alns.append(aln)
  return alns

def _GetMappedResidues(alns):
  """
  :return: Mapping of shared residues in *alns* (with views attached)
           (see :attr:`QSscorer.mapped_residues`).
  :param alns: See :attr:`QSscorer.alignments`
  """
  mapped_residues = dict()
  for aln in alns:
    # prepare dict for c1
    c1 = aln.GetSequence(0).name
    mapped_residues[c1] = dict()
    # get shared residues
    v1, v2 = seq.ViewsFromAlignment(aln)
    for res_1, res_2 in zip(v1.residues, v2.residues):
      r1 = res_1.number.num
      r2 = res_2.number.num
      mapped_residues[c1][r1] = r2

  return mapped_residues

def _GetExtraWeights(contacts, done, mapped_residues):
  """Return sum of extra weights for contacts of chains in set and not in done.
  :return: Tuple (weight_extra_mapped, weight_extra_all).
           weight_extra_mapped only sums if both cX,rX in mapped_residues
           weight_extra_all sums all
  :param contacts: See :func:`GetContacts` for first entity
  :param done: List of (c1, c2, r1, r2) tuples to ignore
  :param mapped_residues: See :func:`_GetMappedResidues`
  """
  weight_extra_mapped = 0
  weight_extra_non_mapped = 0
  for c1 in contacts:
    for c2 in contacts[c1]:
      for r1 in contacts[c1][c2]:
        for r2 in contacts[c1][c2][r1]:
          if (c1, c2, r1, r2) not in done:
            weight = _weight(contacts[c1][c2][r1][r2])
            if     c1 in mapped_residues and r1 in mapped_residues[c1] \
               and c2 in mapped_residues and r2 in mapped_residues[c2]:
              weight_extra_mapped += weight
            else:
              weight_extra_non_mapped += weight
  return weight_extra_mapped, weight_extra_mapped + weight_extra_non_mapped

def _GetScores(contacts_1, contacts_2, mapped_residues, chain_mapping):
  """Get QS scores (see :class:`QSscorer`).

  Note that if some chains are not to be considered at all, they must be removed
  from *contacts_1* / *contacts_2* prior to calling this.

  :param contacts_1: See :func:`GetContacts` for first entity
  :param contacts_2: See :func:`GetContacts` for second entity
  :param mapped_residues: See :func:`_GetMappedResidues`
  :param chain_mapping: Maps any chain name in *mapped_residues* to chain name
                        for *contacts_2*.
  :type chain_mapping:  :class:`dict` (:class:`str` / :class:`str`)
  :return: Tuple (QS_best, QS_global) (see :attr:`QSscorer.best_score` and
           see :attr:`QSscorer.global_score`)
  """
  # keep track of considered contacts as set of (c1,c2,r1,r2) tuples
  done_1 = set()
  done_2 = set()
  weighted_scores = 0
  weight_sum = 0
  # naming cXY, rXY: X = 1/2 for contact, Y = 1/2/T for mapping (T = tmp)
  # -> d1 = contacts_1[c11][c21][r11][r21], d2 = contacts_2[c12][c22][r12][r22]
  for c11 in contacts_1:
    # map to other one
    if c11 not in mapped_residues: continue
    c1T = chain_mapping[c11]
    # second chain
    for c21 in contacts_1[c11]:
      # map to other one
      if c21 not in mapped_residues: continue
      c2T = chain_mapping[c21]
      # flip if needed
      flipped_chains = (c1T > c2T)
      if flipped_chains:
        c12, c22 = c2T, c1T
      else:
        c12, c22 = c1T, c2T
      # skip early if no contacts there
      if c12 not in contacts_2 or c22 not in contacts_2[c12]: continue
      # loop over res.num. in c11
      for r11 in contacts_1[c11][c21]:
        # map to other one and skip if no contacts there
        if r11 not in mapped_residues[c11]: continue
        r1T = mapped_residues[c11][r11]
        # loop over res.num. in c21
        for r21 in contacts_1[c11][c21][r11]:
          # map to other one and skip if no contacts there
          if r21 not in mapped_residues[c21]: continue
          r2T = mapped_residues[c21][r21]
          # flip if needed
          if flipped_chains:
            r12, r22 = r2T, r1T
          else:
            r12, r22 = r1T, r2T
          # skip early if no contacts there
          if r12 not in contacts_2[c12][c22]: continue
          if r22 not in contacts_2[c12][c22][r12]: continue
          # ok now we can finally do the scoring
          d1 = contacts_1[c11][c21][r11][r21]
          d2 = contacts_2[c12][c22][r12][r22]
          weight = _weight(min(d1, d2))
          weighted_scores += weight * (1 - abs(d1 - d2) / 12)
          weight_sum += weight
          # keep track of done ones
          done_1.add((c11, c21, r11, r21))
          done_2.add((c12, c22, r12, r22))

  LogVerbose("Shared contacts: %d" % len(done_1))

  # add the extra weights
  weights_extra_1 = _GetExtraWeights(contacts_1, done_1, mapped_residues)
  mapped_residues_2 = dict()
  for c1 in mapped_residues:
    c2 = chain_mapping[c1]
    mapped_residues_2[c2] = set()
    for r1 in mapped_residues[c1]:
      mapped_residues_2[c2].add(mapped_residues[c1][r1])
  weights_extra_2 = _GetExtraWeights(contacts_2, done_2, mapped_residues_2)
  weight_extra_mapped = weights_extra_1[0] + weights_extra_2[0]
  weight_extra_all    = weights_extra_1[1] + weights_extra_2[1]
  
  # get scores
  denom_best = weight_sum + weight_extra_mapped
  denom_all  = weight_sum + weight_extra_all
  if denom_best == 0:
    QS_best = 0
  else:
    QS_best = weighted_scores / denom_best
  if denom_all == 0:
    QS_global = 0
  else:
    QS_global = weighted_scores / denom_all
  return QS_best, QS_global

def _weight(dist):
  """
  This weight expresses the probability of two residues to interact given the CB
  distance (from Xu et al. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2573399/)
  """
  if dist <= 5.0:
    return 1
  x = (dist-5.0)/4.28;
  return np.exp(-2*x*x)


def _GetQsSuperposition(alns):
  """
  :return: Superposition result based on shared CA atoms in *alns*
           (with views attached) (see :attr:`QSscorer.superposition`).
  :param alns: See :attr:`QSscorer.alignments`
  """
  # check input
  if not alns:
    raise QSscoreError('No successful alignments for superposition!')
  # prepare views
  view_1 = alns[0].GetSequence(0).attached_view.CreateEmptyView()
  view_2 = alns[0].GetSequence(1).attached_view.CreateEmptyView()
  # collect CA from alignments in proper order
  for aln in alns:
    for col in aln:
      res_1 = col.GetResidue(0)
      res_2 = col.GetResidue(1)
      if res_1.IsValid() and res_2.IsValid():
        ca_1 = res_1.FindAtom('CA')
        ca_2 = res_2.FindAtom('CA')
        if ca_1.IsValid() and ca_2.IsValid():
          view_1.AddAtom(ca_1)
          view_2.AddAtom(ca_2)
  # superpose them without chainging entities
  res = mol.alg.SuperposeSVD(view_1, view_2, apply_transform=False)
  return res


def _AddResidue(edi, res, rnum, chain, calpha_only):
  """
  Add residue *res* with res. num. *run* to given *chain* using editor *edi*.
  Either all atoms added or (if *calpha_only*) only CA.
  """
  if calpha_only:
    ca_atom = res.FindAtom('CA')
    if ca_atom.IsValid():
      new_res = edi.AppendResidue(chain, res.name, rnum)
      edi.InsertAtom(new_res, ca_atom.name, ca_atom.pos)
  else:
    new_res = edi.AppendResidue(chain, res.name, rnum)
    for atom in res.atoms:
      edi.InsertAtom(new_res, atom.name, atom.pos)

def _MergeAlignedChains(alns, ent_1, ent_2, calpha_only):
  """
  Create two new entities (based on the alignments attached views) where all
  residues have same numbering (when they're aligned) and they are all pushed to
  a single chain X. Also append extra chains contained in *ent_1* and *ent_2*
  but not contained in *alns*.

  Used for :attr:`QSscorer.lddt_ref` and :attr:`QSscorer.lddt_mdl`

  :param alns: List of alignments with attached views (first sequence: *ent_1*,
               second: *ent_2*). Residue number in single chain is column index
               of current alignment + sum of lengths of all previous alignments
               (order of alignments as in input list).
  :type alns:  See :attr:`QSscorer.alignments`
  :param ent_1: First entity to process.
  :type ent_1:  :class:`~ost.mol.EntityHandle`
  :param ent_2: Second entity to process.
  :type ent_2:  :class:`~ost.mol.EntityHandle`
  :param calpha_only: If True, we only include CA atoms instead of all.
  :type calpha_only:  :class:`bool`

  :return: Tuple of two single chain entities (from *ent_1* and from *ent_2*)
  :rtype:  :class:`tuple` of :class:`~ost.mol.EntityHandle`
  """
  # first new entity
  ent_ren_1 = mol.CreateEntity()
  ed_1 = ent_ren_1.EditXCS()
  new_chain_1 = ed_1.InsertChain('X')
  # second one
  ent_ren_2 = mol.CreateEntity()
  ed_2 = ent_ren_2.EditXCS()
  new_chain_2 = ed_2.InsertChain('X')
  # the alignment already contains sorted chains
  rnum = 0
  chain_done_1 = set()
  chain_done_2 = set()
  for aln in alns:
    chain_done_1.add(aln.GetSequence(0).name)
    chain_done_2.add(aln.GetSequence(1).name)
    for col in aln:
      rnum += 1
      # add valid residues to single chain entities
      res_1 = col.GetResidue(0)
      if res_1.IsValid():
        _AddResidue(ed_1, res_1, rnum, new_chain_1, calpha_only)
      res_2 = col.GetResidue(1)
      if res_2.IsValid():
        _AddResidue(ed_2, res_2, rnum, new_chain_2, calpha_only)
  # extra chains
  for chain in ent_1.chains:
    if chain.name in chain_done_1:
      continue
    for res in chain.residues:
      rnum += 1
      _AddResidue(ed_1, res, rnum, new_chain_1, calpha_only)
  for chain in ent_2.chains:
    if chain.name in chain_done_2:
      continue
    for res in chain.residues:
      rnum += 1
      _AddResidue(ed_2, res, rnum, new_chain_2, calpha_only)
  # get entity names
  ent_ren_1.SetName(aln.GetSequence(0).GetAttachedView().GetName())
  ent_ren_2.SetName(aln.GetSequence(1).GetAttachedView().GetName())
  # connect atoms
  if not conop.GetDefaultLib():
    raise RuntimeError("QSscore computation requires a compound library!")
  pr = conop.RuleBasedProcessor(conop.GetDefaultLib())
  pr.Process(ent_ren_1)  
  ed_1.UpdateICS()
  pr.Process(ent_ren_2)  
  ed_2.UpdateICS()
  return ent_ren_1, ent_ren_2

def _ComputeLDDTScore(ref, mdl):
  """
  :return: lDDT of *mdl* vs *ref* (see :attr:`QSscorer.lddt_score`).
  :param mdl: Reference entity (see :attr:`QSscorer.lddt_mdl`)
  :param ref: Model entity (see :attr:`QSscorer.lddt_ref`)
  """
  # check input
  LogInfo('Reference %s has: %s residues' % (ref.GetName(), ref.residue_count))
  LogInfo('Model %s has: %s residues' % (mdl.GetName(), mdl.residue_count))
  # get lddt score with fixed settings
  lddt_score = mol.alg.LocalDistDiffTest(mdl.Select(''), ref.Select(''),
                                         2., 8., 'lddt')
  LogInfo('lDDT score: %.3f' % lddt_score)
  # add lDDT as B-factor to model
  for r in mdl.residues:
    if r.HasProp('lddt'):
      for a in r.atoms:
        a.SetBFactor(r.GetFloatProp('lddt'))
    else:
      for a in r.atoms:
        a.SetBFactor(0.0)
  
  return lddt_score

# specify public interface
__all__ = ('QSscoreError', 'QSscorer', 'QSscoreEntity', 'FilterContacts',
           'GetContacts')

---