contact_score.py

Go to the documentation of this file.

import itertools
import numpy as np
 
import time
from ost import mol
from ost import geom
from ost import io
 
class ContactEntity:
    """ Helper object for Contact-score computation
    """
    def __init__(self, ent, contact_d = 5.0, contact_mode="aa"):
 
        if contact_mode not in ["aa", "repr"]:
            raise RuntimeError("contact_mode must be in [\"aa\", \"repr\"]")
 
        if contact_mode == "repr":
            for r in ent.residues:
                repr_at = None
                if r.IsPeptideLinking():
                    cb = r.FindAtom("CB")
                    if cb.IsValid():
                        repr_at = cb
                    elif r.GetName() == "GLY":
                        ca = r.FindAtom("CA")
                        if ca.IsValid():
                            repr_at = ca
                elif r.IsNucleotideLinking():
                    c3 = r.FindAtom("C3'")
                    if c3.IsValid():
                        repr_at = c3
                else:
                    raise RuntimeError(f"Only support peptide and nucleotide "
                                       f"residues in \"repr\" contact mode. "
                                       f"Problematic residue: {r}")
                if repr_at is None:
                    raise RuntimeError(f"Residue {r} has no required "
                                       f"representative atom (CB for peptide "
                                       f"residues (CA for GLY) C3' for "
                                       f"nucleotide residues.")
 
        self._contact_mode = contact_mode
 
        if self.contact_modecontact_mode == "aa":
            self._view = ent.CreateFullView()
        elif self.contact_modecontact_mode == "repr":
            pep_query = "(peptide=true and (aname=\"CB\" or (rname=\"GLY\" and aname=\"CA\")))"
            nuc_query = "(nucleotide=True and aname=\"C3'\")"
            self._view = ent.Select(" or ".join([pep_query, nuc_query]))
        self._contact_d = contact_d
 
        # the following attributes will be lazily evaluated
        self._chain_names = None
        self._interacting_chains = None
        self._sequence = dict()
        self._contacts = None
        self._hr_contacts = None
        self._interface_residues = None
        self._hr_interface_residues = None
 
    @property
    def view(self):
        """ The structure depending on *contact_mode*
 
        Full view in case of "aa", view that only contains representative
        atoms in case of "repr".
 
        :type: :class:`ost.mol.EntityView`
        """
        return self._view
 
    @property
    def contact_mode(self):
        """ The contact mode
 
        Can either be "aa", meaning that all atoms are considered to identify
        contacts, or "repr" which only considers distances between
        representative atoms. For peptides thats CB (CA for GLY), for
        nucleotides thats C3'.
 
        :type: :class:`str`
        """
        return self._contact_mode
    
    @property
    def contact_d(self):
        """ Pairwise distance of residues to be considered as contacts
 
        Given at :class:`ContactScorer` construction
 
        :type: :class:`float`
        """
        return self._contact_d
 
    @property
    def chain_names(self):
        """ Chain names in :attr:`~view`
 
        Names are sorted
 
        :type: :class:`list` of :class:`str`
        """
        if self._chain_names is None:
            self._chain_names = sorted([ch.name for ch in self.view.chains])
        return self._chain_names
 
    @property
    def interacting_chains(self):
        """ Pairs of chains in :attr:`~view` with at least one contact
 
        :type: :class:`list` of :class:`tuples`
        """
        if self._interacting_chains is None:
            self._interacting_chains = list(self.contacts.keys())
        return self._interacting_chains
    
    @property
    def contacts(self):
        """ Interchain contacts
 
        Organized as :class:`dict` with key (cname1, cname2) and values being
        a set of tuples with the respective residue indices. 
        cname1 < cname2 evaluates to True.
        """
        if self._contacts is None:
            self._SetupContacts()
        return self._contacts
 
    @property
    def hr_contacts(self):
        """ Human readable interchain contacts
 
        Human readable version of :attr:`~contacts`. Simple list with tuples
        containing two strings specifying the residues in contact. Format:
        <cname>.<rnum>.<ins_code>
        """
        if self._hr_contacts is None:
            self._SetupContacts()
        return self._hr_contacts
 
    @property
    def interface_residues(self):
        """ Interface residues
 
        Residues in each chain that are in contact with any other chain.
        Organized as :class:`dict` with key cname and values the respective
        residue indices in a :class:`set`.
        """
        if self._interface_residues is None:
            self._SetupInterfaceResidues()
        return self._interface_residues
 
    @property
    def hr_interface_residues(self):
        """ Human readable interface residues
 
        Human readable version of :attr:`interface_residues`. :class:`list` of
        strings specifying the interface residues in format:
        <cname>.<rnum>.<ins_code>
        """
        if self._interface_residues is None:
            self._SetupHRInterfaceResidues()
        return self._hr_interface_residues
    
    def GetChain(self, chain_name):
        """ Get chain by name
 
        :param chain_name: Chain in :attr:`~view`
        :type chain_name: :class:`str`
        """ 
        chain = self.view.FindChain(chain_name)
        if not chain.IsValid():
            raise RuntimeError(f"view has no chain named \"{chain_name}\"")
        return chain
 
    def GetSequence(self, chain_name):
        """ Get sequence of chain
 
        Returns sequence of specified chain as raw :class:`str`
 
        :param chain_name: Chain in :attr:`~view`
        :type chain_name: :class:`str`
        """
        if chain_name not in self._sequence:
            ch = self.GetChain(chain_name)
            s = ''.join([r.one_letter_code for r in ch.residues])
            self._sequence[chain_name] = s
        return self._sequence[chain_name]
 
    def _SetupContacts(self):
        self._contacts = dict()
        self._hr_contacts = list()
 
        # set indices relative to full view 
        for ch in self.view.chains:
            for r_idx, r in enumerate(ch.residues):
                r.SetIntProp("contact_idx", r_idx)
 
        residue_lists = list()
        min_res_x = list()
        min_res_y = list()
        min_res_z = list()
        max_res_x = list()
        max_res_y = list()
        max_res_z = list()
        per_res_pos = list()
        min_chain_pos = list()
        max_chain_pos = list()
 
        for cname in self.chain_nameschain_names:
            ch = self.view.FindChain(cname)
            if ch.GetAtomCount() == 0:
                raise RuntimeError(f"Chain without atoms observed: \"{cname}\"")
            residue_lists.append([r for r in ch.residues])
            res_pos = list()
            for r in residue_lists[-1]:
                pos = np.zeros((r.GetAtomCount(), 3))
                for at_idx, at in enumerate(r.atoms):
                    p = at.GetPos()
                    pos[(at_idx, 0)] = p[0]
                    pos[(at_idx, 1)] = p[1]
                    pos[(at_idx, 2)] = p[2]
                res_pos.append(pos)
            min_res_pos = np.vstack([p.min(0) for p in res_pos])
            max_res_pos = np.vstack([p.max(0) for p in res_pos])
            min_res_x.append(min_res_pos[:, 0])
            min_res_y.append(min_res_pos[:, 1])
            min_res_z.append(min_res_pos[:, 2])
            max_res_x.append(max_res_pos[:, 0])
            max_res_y.append(max_res_pos[:, 1])
            max_res_z.append(max_res_pos[:, 2])
            min_chain_pos.append(min_res_pos.min(0))
            max_chain_pos.append(max_res_pos.max(0))
            per_res_pos.append(res_pos)
 
        # operate on squared contact_d (scd) to save some square roots
        scd = self.contact_d * self.contact_d 
 
        for ch1_idx in range(len(self.chain_nameschain_names)):
            for ch2_idx in range(ch1_idx + 1, len(self.chain_nameschain_names)):
                # chains which fulfill the following expressions have no contact
                # within self.contact_d
                if np.max(min_chain_pos[ch1_idx] - max_chain_pos[ch2_idx]) > self.contact_d:
                    continue
                if np.max(min_chain_pos[ch2_idx] - max_chain_pos[ch1_idx]) > self.contact_d:
                    continue
 
                # same thing for residue positions but all at once
                skip_one = np.subtract.outer(min_res_x[ch1_idx], max_res_x[ch2_idx]) > self.contact_d
                skip_one = np.logical_or(skip_one, np.subtract.outer(min_res_y[ch1_idx], max_res_y[ch2_idx]) > self.contact_d)
                skip_one = np.logical_or(skip_one, np.subtract.outer(min_res_z[ch1_idx], max_res_z[ch2_idx]) > self.contact_d)
                skip_two = np.subtract.outer(min_res_x[ch2_idx], max_res_x[ch1_idx]) > self.contact_d
                skip_two = np.logical_or(skip_two, np.subtract.outer(min_res_y[ch2_idx], max_res_y[ch1_idx]) > self.contact_d)
                skip_two = np.logical_or(skip_two, np.subtract.outer(min_res_z[ch2_idx], max_res_z[ch1_idx]) > self.contact_d)
                skip = np.logical_or(skip_one, skip_two.T)
 
                # identify residue pairs for which we cannot exclude a contact
                r1_indices, r2_indices = np.nonzero(np.logical_not(skip))
                ch1_per_res_pos = per_res_pos[ch1_idx]
                ch2_per_res_pos = per_res_pos[ch2_idx]
                for r1_idx, r2_idx in zip(r1_indices, r2_indices):
                    # compute pairwise distances
                    p1 = ch1_per_res_pos[r1_idx]
                    p2 = ch2_per_res_pos[r2_idx]
                    x2 = np.sum(p1**2, axis=1) # (m)
                    y2 = np.sum(p2**2, axis=1) # (n)
                    xy = np.matmul(p1, p2.T) # (m, n)
                    x2 = x2.reshape(-1, 1)
                    squared_distances = x2 - 2*xy + y2 # (m, n)
                    if np.min(squared_distances) <= scd:
                        # its a contact!
                        r1 = residue_lists[ch1_idx][r1_idx]
                        r2 = residue_lists[ch2_idx][r2_idx]
                        cname_key = (self.chain_nameschain_names[ch1_idx], self.chain_nameschain_names[ch2_idx])
                        if cname_key not in self._contacts:
                            self._contacts[cname_key] = set()
                        self._contacts[cname_key].add((r1.GetIntProp("contact_idx"),
                                                       r2.GetIntProp("contact_idx")))
                        rnum1 = r1.GetNumber()
                        hr1 = f"{self.chain_names[ch1_idx]}.{rnum1.num}.{rnum1.ins_code}"
                        rnum2 = r2.GetNumber()
                        hr2 = f"{self.chain_names[ch2_idx]}.{rnum2.num}.{rnum2.ins_code}"
                        self._hr_contacts.append((hr1.strip("\u0000"),
                                                  hr2.strip("\u0000")))
 
 
    def _SetupInterfaceResidues(self):
        self._interface_residues = {cname: set() for cname in self.chain_nameschain_names}
        for k,v in self.contacts.items():
            for item in v:
                self._interface_residues[k[0]].add(item[0])
                self._interface_residues[k[1]].add(item[1])
 
    def _SetupHRInterfaceResidues(self):
        interface_residues = set()
        for item in self.hr_contacts:
            interface_residues.add(item[0])
            interface_residues.add(item[1])
        self._hr_interface_residues = list(interface_residues)
 
 
class ContactScorerResultICS:
    """
    Holds data relevant to compute ics
    """
    def __init__(self, n_trg_contacts, n_mdl_contacts, n_union, n_intersection):
        self._n_trg_contacts = n_trg_contacts
        self._n_mdl_contacts = n_mdl_contacts
        self._n_union = n_union
        self._n_intersection = n_intersection
 
    @property
    def n_trg_contacts(self):
        """ Number of contacts in target
 
        :type: :class:`int`
        """
        return self._n_trg_contacts
 
    @property
    def n_mdl_contacts(self):
        """ Number of contacts in model
 
        :type: :class:`int`
        """
        return self._n_mdl_contacts
 
    @property
    def precision(self):
        """ Precision of model contacts
 
        The fraction of model contacts that are also present in target
 
        :type: :class:`int`
        """
        if self._n_mdl_contacts != 0:
            return self._n_intersection / self._n_mdl_contacts
        else:
            return 0.0
 
    @property
    def recall(self):
        """ Recall of model contacts
 
        The fraction of target contacts that are also present in model
 
        :type: :class:`int`
        """
        if self._n_trg_contacts != 0:
            return self._n_intersection / self._n_trg_contacts
        else:
            return 0.0
 
    @property
    def ics(self):
        """ The Interface Contact Similarity score (ICS)
 
        Combination of :attr:`precision` and :attr:`recall` using the F1-measure
 
        :type: :class:`float`
        """
        p = self.precision
        r = self.recall
        nominator = p*r
        denominator = p + r
        if denominator != 0.0:
            return 2*nominator/denominator
        else:
            return 0.0
 
class ContactScorerResultIPS:
    """
    Holds data relevant to compute ips
    """
    def __init__(self, n_trg_int_res, n_mdl_int_res, n_union, n_intersection):
        self._n_trg_int_res = n_trg_int_res
        self._n_mdl_int_res = n_mdl_int_res
        self._n_union = n_union
        self._n_intersection = n_intersection
 
    @property
    def n_trg_int_res(self):
        """ Number of interface residues in target
 
        :type: :class:`int`
        """
        return self._n_trg_int_res
 
    @property
    def n_mdl_int_res(self):
        """ Number of interface residues in model
 
        :type: :class:`int`
        """
        return self._n_mdl_int_res
 
    @property
    def precision(self):
        """ Precision of model interface residues
 
        The fraction of model interface residues that are also interface
        residues in target
 
        :type: :class:`int`
        """
        if self._n_mdl_int_res != 0:
            return self._n_intersection / self._n_mdl_int_res
        else:
            return 0.0
 
    @property
    def recall(self):
        """ Recall of model interface residues
 
        The fraction of target interface residues that are also interface
        residues in model
 
        :type: :class:`int`
        """
        if self._n_trg_int_res != 0:
            return self._n_intersection / self._n_trg_int_res
        else:
            return 0.0
 
    @property
    def ips(self):
        """ The Interface Patch Similarity score (IPS)
 
        Jaccard coefficient of interface residues in model/target.
        Technically thats :attr:`intersection`/:attr:`union` 
 
        :type: :class:`float`
        """
        if(self._n_union > 0):
            return self._n_intersection/self._n_union
        return 0.0
 
class ContactScorer:
    """ Helper object to compute Contact scores
 
    Tightly integrated into the mechanisms from the chain_mapping module.
    The prefered way to derive an object of type :class:`ContactScorer` is
    through the static constructor: :func:`~FromMappingResult`.
 
    Usage is the same as for :class:`ost.mol.alg.QSScorer`
    """
 
    def __init__(self, target, chem_groups, model, alns,
                 contact_mode="aa", contact_d=5.0):
        self._cent1 = ContactEntity(target, contact_mode = contact_mode,
                                    contact_d = contact_d)
        # ensure that target chain names match the ones in chem_groups
        chem_group_ch_names = list(itertools.chain.from_iterable(chem_groups))
        if self._cent1.chain_names != sorted(chem_group_ch_names):
            raise RuntimeError(f"Expect exact same chain names in chem_groups "
                               f"and in target (which is processed to only "
                               f"contain peptides/nucleotides). target: "
                               f"{self._cent1.chain_names}, chem_groups: "
                               f"{chem_group_ch_names}")
 
        self._chem_groups = chem_groups
        self._cent2 = ContactEntity(model, contact_mode = contact_mode,
                                    contact_d = contact_d)
        self._alns = alns
 
        # cache for mapped interface scores
        # key: tuple of tuple ((qsent1_ch1, qsent1_ch2),
        #                     ((qsent2_ch1, qsent2_ch2))
        # value: tuple with four numbers required for computation of
        #        per-interface scores.
        #        The first two are relevant for ICS, the others for per
        #        interface IPS.
        #        1: n_union_contacts
        #        2: n_intersection_contacts
        #        3: n_union_interface_residues
        #        4: n_intersection_interface_residues
        self._mapped_cache_interface = dict()
 
        # cache for mapped single chain scores
        # for interface residues of single chains
        # key: tuple: (qsent1_ch, qsent2_ch)
        # value: tuple with two numbers required for computation of IPS
        #        1: n_union
        #        2: n_intersection
        self._mapped_cache_sc = dict()
 
    @staticmethod
    def FromMappingResult(mapping_result, contact_mode="aa", contact_d = 5.0):
        """ The preferred way to get a :class:`ContactScorer`
 
        Static constructor that derives an object of type :class:`ContactScorer`
        using a :class:`ost.mol.alg.chain_mapping.MappingResult`
 
        :param mapping_result: Data source
        :type mapping_result: :class:`ost.mol.alg.chain_mapping.MappingResult`
        """
        contact_scorer = ContactScorer(mapping_result.target,
                                       mapping_result.chem_groups,
                                       mapping_result.model,
                                       mapping_result.alns,
                                       contact_mode = contact_mode,
                                       contact_d = contact_d)
        return contact_scorer
 
    @property
    def cent1(self):
        """ Represents *target*
 
        :type: :class:`ContactEntity`
        """
        return self._cent1
 
    @property
    def chem_groups(self):
        """ Groups of chemically equivalent chains in *target*
 
        Provided at object construction
 
        :type: :class:`list` of :class:`list` of :class:`str`
        """
        return self._chem_groups
 
    @property
    def cent2(self):
        """ Represents *model*
 
        :type: :class:`ContactEntity`
        """
        return self._cent2
 
    @property
    def alns(self):
        """ Alignments between chains in :attr:`~cent1` and :attr:`~cent2`
 
        Provided at object construction. Each alignment is accessible with
        ``alns[(t_chain,m_chain)]``. First sequence is the sequence of the
        respective chain in :attr:`~cent1`, second sequence the one from
        :attr:`~cent2`.
 
        :type: :class:`dict` with key: :class:`tuple` of :class:`str`, value:
               :class:`ost.seq.AlignmentHandle`
        """
        return self._alns
 
    def ScoreICS(self, mapping, check=True):
        """ Computes ICS given chain mapping
 
        Again, the preferred way is to get *mapping* is from an object
        of type :class:`ost.mol.alg.chain_mapping.MappingResult`.
 
        :param mapping: see 
                        :attr:`ost.mol.alg.chain_mapping.MappingResult.mapping`
        :type mapping: :class:`list` of :class:`list` of :class:`str`
        :param check: Perform input checks, can be disabled for speed purposes
                      if you know what you're doing.
        :type check: :class:`bool`
        :returns: Result object of type :class:`ContactScorerResultICS`
        """
 
        if check:
            # ensure that dimensionality of mapping matches self.chem_groups
            if len(self.chem_groupschem_groups) != len(mapping):
                raise RuntimeError("Dimensions of self.chem_groups and mapping "
                                   "must match")
            for a,b in zip(self.chem_groupschem_groups, mapping):
                if len(a) != len(b):
                    raise RuntimeError("Dimensions of self.chem_groups and "
                                       "mapping must match")
            # ensure that chain names in mapping are all present in cent2
            for name in itertools.chain.from_iterable(mapping):
                if name is not None and name not in self.cent2.chain_names:
                    raise RuntimeError(f"Each chain in mapping must be present "
                                       f"in self.cent2. No match for "
                                       f"\"{name}\"")
 
        flat_mapping = dict()
        for a, b in zip(self.chem_groupschem_groups, mapping):
            flat_mapping.update({x: y for x, y in zip(a, b) if y is not None})
 
        return self.ICSFromFlatMapping(flat_mapping)
 
    def ScoreICSInterface(self, trg_ch1, trg_ch2, mdl_ch1, mdl_ch2):
        """ Computes ICS scores only considering one interface
 
        This only works for interfaces that are computed in :func:`Score`, i.e.
        interfaces for which the alignments are set up correctly.
 
        :param trg_ch1: Name of first interface chain in target
        :type trg_ch1: :class:`str`
        :param trg_ch2: Name of second interface chain in target
        :type trg_ch2: :class:`str`
        :param mdl_ch1: Name of first interface chain in model
        :type mdl_ch1: :class:`str`
        :param mdl_ch2: Name of second interface chain in model
        :type mdl_ch2: :class:`str`
        :returns: Result object of type :class:`ContactScorerResultICS`
        :raises: :class:`RuntimeError` if no aln for trg_ch1/mdl_ch1 or
                 trg_ch2/mdl_ch2 is available.
        """
        if (trg_ch1, mdl_ch1) not in self.alns:
            raise RuntimeError(f"No aln between trg_ch1 ({trg_ch1}) and "
                               f"mdl_ch1 ({mdl_ch1}) available. Did you "
                               f"construct the QSScorer object from a "
                               f"MappingResult and are trg_ch1 and mdl_ch1 "
                               f"mapped to each other?")
        if (trg_ch2, mdl_ch2) not in self.alns:
            raise RuntimeError(f"No aln between trg_ch1 ({trg_ch1}) and "
                               f"mdl_ch1 ({mdl_ch1}) available. Did you "
                               f"construct the QSScorer object from a "
                               f"MappingResult and are trg_ch1 and mdl_ch1 "
                               f"mapped to each other?")
        trg_int = (trg_ch1, trg_ch2)
        mdl_int = (mdl_ch1, mdl_ch2)
        trg_int_r = (trg_ch2, trg_ch1)
        mdl_int_r = (mdl_ch2, mdl_ch1)
 
        if trg_int in self.cent1.contacts:
            n_trg = len(self.cent1.contacts[trg_int])
        elif trg_int_r in self.cent1.contacts:
            n_trg = len(self.cent1.contacts[trg_int_r])
        else:
            n_trg = 0
 
        if mdl_int in self.cent2.contacts:
            n_mdl = len(self.cent2.contacts[mdl_int])
        elif mdl_int_r in self.cent2.contacts:
            n_mdl = len(self.cent2.contacts[mdl_int_r])
        else:
            n_mdl = 0
 
        n_union, n_intersection, _, _ = self._MappedInterfaceScores(trg_int, mdl_int)
        return ContactScorerResultICS(n_trg, n_mdl, n_union, n_intersection)
 
    def ICSFromFlatMapping(self, flat_mapping):
        """ Same as :func:`ScoreICS` but with flat mapping
 
        :param flat_mapping: Dictionary with target chain names as keys and
                             the mapped model chain names as value
        :type flat_mapping: :class:`dict` with :class:`str` as key and value
        :returns: Result object of type :class:`ContactScorerResultICS`
        """
        n_trg = sum([len(x) for x in self.cent1.contacts.values()])
        n_mdl = sum([len(x) for x in self.cent2.contacts.values()])
        n_union = 0
        n_intersection = 0
 
        processed_cent2_interfaces = set()
        for int1 in self.cent1.interacting_chains:
            if int1[0] in flat_mapping and int1[1] in flat_mapping:
                int2 = (flat_mapping[int1[0]], flat_mapping[int1[1]])
                a, b, _, _ = self._MappedInterfaceScores(int1, int2)
                n_union += a
                n_intersection += b
                processed_cent2_interfaces.add((min(int2), max(int2)))
 
        # process interfaces that only exist in qsent2
        r_flat_mapping = {v:k for k,v in flat_mapping.items()} # reverse mapping
        for int2 in self.cent2.interacting_chains:
            if int2 not in processed_cent2_interfaces:
                if int2[0] in r_flat_mapping and int2[1] in r_flat_mapping:
                    int1 = (r_flat_mapping[int2[0]], r_flat_mapping[int2[1]])
                    a, b, _, _ = self._MappedInterfaceScores(int1, int2)
                    n_union += a
                    n_intersection += b
 
        return ContactScorerResultICS(n_trg, n_mdl,
                                      n_union, n_intersection)
 
    def ScoreIPS(self, mapping, check=True):
        """ Computes IPS given chain mapping
 
        Again, the preferred way is to get *mapping* is from an object
        of type :class:`ost.mol.alg.chain_mapping.MappingResult`.
 
        :param mapping: see 
                        :attr:`ost.mol.alg.chain_mapping.MappingResult.mapping`
        :type mapping: :class:`list` of :class:`list` of :class:`str`
        :param check: Perform input checks, can be disabled for speed purposes
                      if you know what you're doing.
        :type check: :class:`bool`
        :returns: Result object of type :class:`ContactScorerResultIPS`
        """
 
        if check:
            # ensure that dimensionality of mapping matches self.chem_groups
            if len(self.chem_groupschem_groups) != len(mapping):
                raise RuntimeError("Dimensions of self.chem_groups and mapping "
                                   "must match")
            for a,b in zip(self.chem_groupschem_groups, mapping):
                if len(a) != len(b):
                    raise RuntimeError("Dimensions of self.chem_groups and "
                                       "mapping must match")
            # ensure that chain names in mapping are all present in cent2
            for name in itertools.chain.from_iterable(mapping):
                if name is not None and name not in self.cent2.chain_names:
                    raise RuntimeError(f"Each chain in mapping must be present "
                                       f"in self.cent2. No match for "
                                       f"\"{name}\"")
 
        flat_mapping = dict()
        for a, b in zip(self.chem_groupschem_groups, mapping):
            flat_mapping.update({x: y for x, y in zip(a, b) if y is not None})
 
        return self.IPSFromFlatMapping(flat_mapping)
 
    def ScoreIPSInterface(self, trg_ch1, trg_ch2, mdl_ch1, mdl_ch2):
        """ Computes IPS scores only considering one interface
 
        This only works for interfaces that are computed in :func:`Score`, i.e.
        interfaces for which the alignments are set up correctly.
 
        :param trg_ch1: Name of first interface chain in target
        :type trg_ch1: :class:`str`
        :param trg_ch2: Name of second interface chain in target
        :type trg_ch2: :class:`str`
        :param mdl_ch1: Name of first interface chain in model
        :type mdl_ch1: :class:`str`
        :param mdl_ch2: Name of second interface chain in model
        :type mdl_ch2: :class:`str`
        :returns: Result object of type :class:`ContactScorerResultIPS`
        :raises: :class:`RuntimeError` if no aln for trg_ch1/mdl_ch1 or
                 trg_ch2/mdl_ch2 is available.
        """
        if (trg_ch1, mdl_ch1) not in self.alns:
            raise RuntimeError(f"No aln between trg_ch1 ({trg_ch1}) and "
                               f"mdl_ch1 ({mdl_ch1}) available. Did you "
                               f"construct the QSScorer object from a "
                               f"MappingResult and are trg_ch1 and mdl_ch1 "
                               f"mapped to each other?")
        if (trg_ch2, mdl_ch2) not in self.alns:
            raise RuntimeError(f"No aln between trg_ch1 ({trg_ch1}) and "
                               f"mdl_ch1 ({mdl_ch1}) available. Did you "
                               f"construct the QSScorer object from a "
                               f"MappingResult and are trg_ch1 and mdl_ch1 "
                               f"mapped to each other?")
        trg_int = (trg_ch1, trg_ch2)
        mdl_int = (mdl_ch1, mdl_ch2)
        trg_int_r = (trg_ch2, trg_ch1)
        mdl_int_r = (mdl_ch2, mdl_ch1)
 
        trg_contacts = None
        if trg_int in self.cent1.contacts:
            trg_contacts = self.cent1.contacts[trg_int]
        elif trg_int_r in self.cent1.contacts:
            trg_contacts = self.cent1.contacts[trg_int_r]
 
        if trg_contacts is None:
            n_trg = 0
        else:
            n_trg = len(set([x[0] for x in trg_contacts]))
            n_trg += len(set([x[1] for x in trg_contacts]))
 
        mdl_contacts = None
        if mdl_int in self.cent2.contacts:
            mdl_contacts = self.cent2.contacts[mdl_int]
        elif mdl_int_r in self.cent2.contacts:
            mdl_contacts = self.cent2.contacts[mdl_int_r]
 
        if mdl_contacts is None:
            n_mdl = 0
        else:
            n_mdl = len(set([x[0] for x in mdl_contacts]))
            n_mdl += len(set([x[1] for x in mdl_contacts]))
 
        _, _, n_union, n_intersection = self._MappedInterfaceScores(trg_int, mdl_int)
        return ContactScorerResultIPS(n_trg, n_mdl, n_union, n_intersection)
 
 
    def IPSFromFlatMapping(self, flat_mapping):
        """ Same as :func:`ScoreIPS` but with flat mapping
 
        :param flat_mapping: Dictionary with target chain names as keys and
                             the mapped model chain names as value
        :type flat_mapping: :class:`dict` with :class:`str` as key and value
        :returns: Result object of type :class:`ContactScorerResultIPS`
        """
        n_trg = sum([len(x) for x in self.cent1.interface_residues.values()])
        n_mdl = sum([len(x) for x in self.cent2.interface_residues.values()])
        n_union = 0
        n_intersection = 0
 
        processed_cent2_chains = set()
        for trg_ch in self.cent1.chain_names:
            if trg_ch in flat_mapping:
                a, b = self._MappedSCScores(trg_ch, flat_mapping[trg_ch])
                n_union += a
                n_intersection += b
                processed_cent2_chains.add(flat_mapping[trg_ch])
            else:
                n_union += len(self.cent1.interface_residues[trg_ch])
 
        for mdl_ch in self._cent2.chain_names:
            if mdl_ch not in processed_cent2_chains:
                n_union += len(self.cent2.interface_residues[mdl_ch])
 
        return ContactScorerResultIPS(n_trg, n_mdl,
                                      n_union, n_intersection)
 
 
    def _MappedInterfaceScores(self, int1, int2):
        key_one = (int1, int2)
        if key_one in self._mapped_cache_interface:
            return self._mapped_cache_interface[key_one]
        key_two = ((int1[1], int1[0]), (int2[1], int2[0]))
        if key_two in self._mapped_cache_interface:
            return self._mapped_cache_interface[key_two]
 
        a, b, c, d = self._InterfaceScores(int1, int2)
        self._mapped_cache_interface[key_one] = (a, b, c, d)
        return (a, b, c, d)
 
    def _InterfaceScores(self, int1, int2):
        if int1 in self.cent1.contacts:
            ref_contacts = self.cent1.contacts[int1]
        elif (int1[1], int1[0]) in self.cent1.contacts:
            ref_contacts = self.cent1.contacts[(int1[1], int1[0])]
            # need to reverse contacts
            ref_contacts = set([(x[1], x[0]) for x in ref_contacts])
        else:
            ref_contacts = set() # no contacts at all
 
        if int2 in self.cent2.contacts:
            mdl_contacts = self.cent2.contacts[int2]
        elif (int2[1], int2[0]) in self.cent2.contacts:
            mdl_contacts = self.cent2.contacts[(int2[1], int2[0])]
            # need to reverse contacts
            mdl_contacts = set([(x[1], x[0]) for x in mdl_contacts])
        else:
            mdl_contacts = set() # no contacts at all
 
        # indices in contacts lists are specific to the respective
        # structures, need manual mapping from alignments
        ch1_aln = self.alns[(int1[0], int2[0])]
        ch2_aln = self.alns[(int1[1], int2[1])]
        mapped_ref_contacts = set()
        mapped_mdl_contacts = set()
        for c in ref_contacts:
            mapped_c = (ch1_aln.GetPos(0, c[0]), ch2_aln.GetPos(0, c[1]))
            mapped_ref_contacts.add(mapped_c)
        for c in mdl_contacts:
            mapped_c = (ch1_aln.GetPos(1, c[0]), ch2_aln.GetPos(1, c[1]))
            mapped_mdl_contacts.add(mapped_c)
 
        contact_union = len(mapped_ref_contacts.union(mapped_mdl_contacts))
        contact_intersection = len(mapped_ref_contacts.intersection(mapped_mdl_contacts))
 
        # above, we computed the union and intersection on actual
        # contacts. Here, we do the same on interface residues
 
        # process interface residues of chain one in interface
        tmp_ref = set([x[0] for x in mapped_ref_contacts])
        tmp_mdl = set([x[0] for x in mapped_mdl_contacts])
        intres_union = len(tmp_ref.union(tmp_mdl))
        intres_intersection = len(tmp_ref.intersection(tmp_mdl))
 
        # process interface residues of chain two in interface
        tmp_ref = set([x[1] for x in mapped_ref_contacts])
        tmp_mdl = set([x[1] for x in mapped_mdl_contacts])
        intres_union += len(tmp_ref.union(tmp_mdl))
        intres_intersection += len(tmp_ref.intersection(tmp_mdl))
 
        return (contact_union, contact_intersection,
                intres_union, intres_intersection)
 
    def _MappedSCScores(self, ref_ch, mdl_ch):
        if (ref_ch, mdl_ch) in self._mapped_cache_sc:
            return self._mapped_cache_sc[(ref_ch, mdl_ch)]
        n_union, n_intersection = self._SCScores(ref_ch, mdl_ch)
        self._mapped_cache_sc[(ref_ch, mdl_ch)] = (n_union, n_intersection)
        return (n_union, n_intersection)
 
    def _SCScores(self, ch1, ch2):
        ref_int_res = self.cent1.interface_residues[ch1]
        mdl_int_res = self.cent2.interface_residues[ch2]
        aln = self.alns[(ch1, ch2)]
        mapped_ref_int_res = set()
        mapped_mdl_int_res = set()
        for r_idx in ref_int_res:
            mapped_ref_int_res.add(aln.GetPos(0, r_idx))
        for r_idx in mdl_int_res:
            mapped_mdl_int_res.add(aln.GetPos(1, r_idx))
        return(len(mapped_ref_int_res.union(mapped_mdl_int_res)),
               len(mapped_ref_int_res.intersection(mapped_mdl_int_res)))
 
# specify public interface
__all__ = ('ContactEntity', 'ContactScorerResultICS', 'ContactScorerResultIPS', 'ContactScorer')