chemosim-lab / ProLIF

Interaction Fingerprints for protein-ligand complexes and more
https://prolif.readthedocs.io
Apache License 2.0
361 stars 68 forks source link

VdWContact with hydrophobic SMARTS #179

Closed Dan-Burns closed 9 months ago

Dan-Burns commented 9 months ago

Great work on Prolif!

I'd like to apply the tolerance option for hydrophobic interactions. Can the "detect" method in VdWContact be limited to the atoms picked up in the Hydrophobic classes SMARTS string?

Thank you, Dan

cbouy commented 9 months ago

Not directly, you'd have to pretty much copy-paste both implementations to have the VdW threshold applied to the substructure matches of the Hydrophobic SMARTS, something like this:

from itertools import product

from prolif.interactions import Hydrophobic, VdWContact

class HydrophobicVdW(Hydrophobic, VdWContact):
    def __init__(self, tolerance=0):
        Hydrophobic.__init__(self)
        VdWContact.__init__(self, tolerance=tolerance)

    def detect(self, lig_res, prot_res):
        # match hydrophobic SMARTS
        lig_matches = lig_res.GetSubstructMatches(self.lig_pattern)
        prot_matches = prot_res.GetSubstructMatches(self.prot_pattern)
        if lig_matches and prot_matches:
            # satisfy VdW radii sum cutoff
            lxyz = lig_res.GetConformer()
            rxyz = prot_res.GetConformer()
            for lig_match, prot_match in product(lig_matches, prot_matches):
                la = lig_res.GetAtomWithIdx(lig_match[0])
                ra = prot_res.GetAtomWithIdx(prot_match[0])
                lig = la.GetSymbol()
                res = ra.GetSymbol()
                elements = frozenset((lig, res))
                try:
                    vdw = self._vdw_cache[elements]
                except KeyError:
                    vdw = self.vdwradii[lig] + self.vdwradii[res] + self.tolerance
                    self._vdw_cache[elements] = vdw
                dist = lxyz.GetAtomPosition(lig_match[0]).Distance(
                    rxyz.GetAtomPosition(prot_match[0])
                )
                if dist <= vdw:
                    yield self.metadata(
                        lig_res, prot_res, lig_match, prot_match, distance=dist
                    )

You might have to put this in a separate python file, and then import it and use like

from custom import HydrophobicVdW

fp = plf.Fingerprint(["HydrophobicVdW"])
fp.run(u.trajectory[::10], ligand_selection, protein_selection)

Feel free to reopen the issue if there are other problems/questions