[med-svn] [Git][med-team/plip][master] 6 commits: New upstream version 2.4.0+dfsg

Alexandre Mestiashvili (@mestia) gitlab at salsa.debian.org
Tue Jan 14 14:07:48 GMT 2025



Alexandre Mestiashvili pushed to branch master at Debian Med / plip


Commits:
61a58349 by Alexandre Mestiashvili at 2025-01-10T14:29:38+01:00
New upstream version 2.4.0+dfsg
- - - - -
bec9fb6d by Alexandre Mestiashvili at 2025-01-10T14:29:38+01:00
Update upstream source from tag 'upstream/2.4.0+dfsg'

Update to upstream version '2.4.0+dfsg'
with Debian dir 933dcfd2ff88b579a7a7ffedefb57df1f43945d1
- - - - -
1e768ecf by Alexandre Mestiashvili at 2025-01-10T14:30:10+01:00
Update changelog

- - - - -
ccd9b03c by Alexandre Mestiashvili at 2025-01-14T14:48:20+01:00
Add patch removing custom build option for broken openbabel pip bindings

See also: https://github.com/openbabel/openbabel/issues/2408

- - - - -
c087e0d9 by Alexandre Mestiashvili at 2025-01-14T14:57:24+01:00
d/control, bump Standards-Version to 4.7.0

- - - - -
987142df by Alexandre Mestiashvili at 2025-01-14T15:00:17+01:00
Update changelog, prepare for release

- - - - -


14 changed files:

- .gitignore
- CHANGES.txt
- debian/changelog
- debian/control
- debian/patches/series
- + debian/patches/setup.patch
- plip/basic/config.py
- plip/basic/supplemental.py
- plip/plipcmd.py
- plip/structure/detection.py
- plip/structure/preparation.py
- plip/test/test_literature_validated.py
- plip/visualization/visualize.py
- setup.py


Changes:

=====================================
.gitignore
=====================================
@@ -48,3 +48,4 @@ target/
 # Other
 *~
 .nfs*
+tests/
\ No newline at end of file


=====================================
CHANGES.txt
=====================================
@@ -1,5 +1,16 @@
 Changelog
 ---------
+# 2.4.0
+* new "--chains" flag to enable detection of interactions between protein chains by @PhiCMS and @snbolz
+* update setup.py, attempt to fix and install broken python openbabel bindings 3.1.1.1
+* update Dockerfile
+
+# 2.3.1
+* fixes in interaction detection with peptide ligands by @kalinni in #153
+* fix in hydrogen bond acceptor identification by @kalinni in #151
+* fixes in waterbridge detection and filtering by @kalinni in #152
+* fix NumPy DeprecationWarning when using numpy.uint32 with negative numbers by @QY0831 in #150
+
 # 2.3.0
 * fixes an issue that caused encoding errors to be thrown
 * salt bridge interaction improved


=====================================
debian/changelog
=====================================
@@ -1,3 +1,11 @@
+plip (2.4.0+dfsg-1) unstable; urgency=medium
+
+  * New upstream version 2.4.0+dfsg
+  * Add patch removing custom build option for broken openbabel pip bindings
+  * d/control, bump Standards-Version to 4.7.0
+
+ -- Alexandre Mestiashvili <mestia at debian.org>  Tue, 14 Jan 2025 14:58:59 +0100
+
 plip (2.3.0+dfsg-2) unstable; urgency=medium
 
   * Team upload.


=====================================
debian/control
=====================================
@@ -11,7 +11,7 @@ Build-Depends: debhelper-compat (= 13),
                python3-numpy,
                python3-openbabel,
                python3-setuptools
-Standards-Version: 4.6.2
+Standards-Version: 4.7.0
 Vcs-Browser: https://salsa.debian.org/med-team/plip
 Vcs-Git: https://salsa.debian.org/med-team/plip.git
 Homepage: https://projects.biotec.tu-dresden.de/plip-web/plip/


=====================================
debian/patches/series
=====================================
@@ -1 +1,2 @@
 revert_imagemagick_changes.patch
+setup.patch


=====================================
debian/patches/setup.patch
=====================================
@@ -0,0 +1,16 @@
+Description: Do not call custom cmdclass for building openbabel bindings
+Author: Alex Mestiashvili <mestia at debian.org>
+--- plip.orig/setup.py
++++ plip/setup.py
+@@ -90,10 +90,10 @@
+     license='GPLv2',
+     packages=find_packages(),
+     scripts=['plip/plipcmd.py'],
+-    cmdclass={'build': CustomBuild, 'build_ext': CustomBuildExt, 'install': CustomInstall},
+     install_requires=[
+         'numpy',
+         'lxml'
++        'openbabel'
+         ],
+     entry_points={
+         "console_scripts": [


=====================================
plip/basic/config.py
=====================================
@@ -1,4 +1,4 @@
-__version__ = '2.3.0'
+__version__ = '2.4.0'
 __maintainer__ = 'PharmAI GmbH (2020-2021) - www.pharm.ai - hello at pharm.ai'
 __citation_information__ = "Adasme,M. et al. PLIP 2021: expanding the scope of the protein-ligand interaction profiler to DNA and RNA. " \
                            "Nucl. Acids Res. (05 May 2021), gkab294. doi: 10.1093/nar/gkab294"
@@ -25,12 +25,15 @@ NOFIX = False  # Turn off fixing of errors in PDB files
 NOFIXFILE = False  # Turn off writing to files for fixed PDB structures
 PEPTIDES = []  # Definition which chains should be considered as peptide ligands
 INTRA = None
+RESIDUES = {}
 KEEPMOD = False
 DNARECEPTOR = False
 OUTPUTFILENAME = "report"  # Naming for the TXT and XML report files
 NOPDBCANMAP = False  # Skip calculation of mapping canonical atom order: PDB atom order
 NOHYDRO = False  # Do not add hydrogen bonds (in case already present in the structure)
 MODEL = 1  # The model to be selected for multi-model structures (default = 1).
+CHAINS = None # Define chains for protein-protein interaction detection
+
 
 # Configuration file for Protein-Ligand Interaction Profiler (PLIP)
 # Set thresholds for detection of interactions


=====================================
plip/basic/supplemental.py
=====================================
@@ -55,6 +55,18 @@ def whichchain(atom):
     return atom.GetResidue().GetChain() if atom.GetResidue() is not None else None
 
 
+def residue_belongs_to_receptor(res, config):
+    """tests whether the residue is defined as receptor and is not part of a peptide or residue ligand."""
+    if config.CHAINS:
+        if config.CHAINS[0] and config.CHAINS[1]:  # if receptor and ligand chains were given
+            return res.GetChain() in config.CHAINS[0] and res.GetChain() not in config.CHAINS[1]
+            # True if residue is part of receptor chains and not of ligand chains
+        if config.CHAINS[1]:  # if only ligand chains were given
+            return res.GetChain() not in config.CHAINS[1]  # True if residue is not part of ligand chains
+        return False  # if only receptor chains were given or both is empty
+    return res.GetChain() not in config.PEPTIDES  # True if residue is not part of peptide ligand.
+
+
 #########################
 # Mathematical operations
 #########################
@@ -358,7 +370,7 @@ def int32_to_negative(int32):
     if int32 == 4294967295:  # Special case in some structures (note, this is just a workaround)
         return -1
     for i in range(-1000, -1):
-        dct[np.uint32(i)] = i
+        dct[int(np.array(i).astype(np.uint32))] = i
     if int32 in dct:
         return dct[int32]
     else:


=====================================
plip/plipcmd.py
=====================================
@@ -10,6 +10,7 @@ import logging
 import multiprocessing
 import os
 import sys
+import ast
 from argparse import ArgumentParser
 from collections import namedtuple
 
@@ -38,9 +39,19 @@ def threshold_limiter(aparser, arg):
         aparser.error("All thresholds have to be values larger than zero.")
     return arg
 
+def residue_list(input_string):
+    """Parse mix of residue numbers and ranges passed with the --residues flag into one list"""
+    result = []
+    for part in input_string.split(','):
+        if '-' in part:
+            start, end = map(int, part.split('-'))
+            result.extend(range(start, end + 1))
+        else:
+            result.append(int(part))
+    return result
 
 def process_pdb(pdbfile, outpath, as_string=False, outputprefix='report'):
-    """Analysis of a single PDB file. Can generate textual reports XML, PyMOL session files and images as output."""
+    """Analysis of a single PDB file with optional chain filtering."""
     if not as_string:
         pdb_file_name = pdbfile.split('/')[-1]
         startmessage = f'starting analysis of {pdb_file_name}'
@@ -50,7 +61,6 @@ def process_pdb(pdbfile, outpath, as_string=False, outputprefix='report'):
     mol = PDBComplex()
     mol.output_path = outpath
     mol.load_pdb(pdbfile, as_string=as_string)
-    # @todo Offers possibility for filter function from command line (by ligand chain, position, hetid)
     for ligand in mol.ligands:
         mol.characterize_complex(ligand)
 
@@ -116,7 +126,7 @@ def remove_duplicates(slist):
     return unique
 
 
-def run_analysis(inputstructs, inputpdbids):
+def run_analysis(inputstructs, inputpdbids, chains=None):
     """Main function. Calls functions for processing, report generation and visualization."""
     pdbid, pdbpath = None, None
     # @todo For multiprocessing, implement better stacktracing for errors
@@ -127,16 +137,16 @@ def run_analysis(inputstructs, inputpdbids):
     output_prefix = config.OUTPUTFILENAME
 
     if inputstructs is not None:  # Process PDB file(s)
-        num_structures = len(inputstructs)
+        num_structures = len(inputstructs)  # @question: how can it become more than one file? The tilde_expansion function does not consider this case.
         inputstructs = remove_duplicates(inputstructs)
         read_from_stdin = False
         for inputstruct in inputstructs:
-            if inputstruct == '-':
+            if inputstruct == '-':  # @expl: when user gives '-' as input, pdb file is read from stdin
                 inputstruct = sys.stdin.read()
                 read_from_stdin = True
                 if config.RAWSTRING:
                     if sys.version_info < (3,):
-                        inputstruct = bytes(inputstruct).decode('unicode_escape')
+                        inputstruct = bytes(inputstruct).decode('unicode_escape')  # @expl: in Python2, the bytes object is just a string.
                     else:
                         inputstruct = bytes(inputstruct, 'utf8').decode('unicode_escape')
             else:
@@ -218,6 +228,9 @@ def main():
                             help="Allows to define one or multiple chains as peptide ligands or to detect inter-chain contacts",
                             nargs="+")
     ligandtype.add_argument("--intra", dest="intra", help="Allows to define one chain to analyze intra-chain contacts.")
+    parser.add_argument("--residues", dest="residues", default=[], nargs="+",
+                        help="""Allows to specify which residues of the chain(s) should be considered as peptide ligands.
+                        Give single residues (separated with comma) or ranges (with dash) or both, for several chains separate selections with one space""")
     parser.add_argument("--keepmod", dest="keepmod", default=False,
                         help="Keep modified residues as ligands",
                         action="store_true")
@@ -241,7 +254,19 @@ def main():
     for t in thresholds:
         parser.add_argument('--%s' % t.name, dest=t.name, type=lambda val: threshold_limiter(parser, val),
                             help=argparse.SUPPRESS)
+
+    # Add argument to define receptor and ligand chains
+    parser.add_argument("--chains", dest="chains", type=str,
+                        help="Specify chains as receptor/ligand groups, e.g., '[['A'], ['B']]'. "
+                             "Use format [['A'], ['B', 'C']] to define A as receptor, and B, C as ligands.")
+
+
     arguments = parser.parse_args()
+    # make sure, residues is only used together with --inter (could be expanded to --intra in the future)
+    if arguments.residues and not (arguments.peptides or arguments.intra):
+        parser.error("The --residues option requires specification of a chain with --inter or --peptide")
+    if arguments.residues and len(arguments.residues)!=len(arguments.peptides):
+        parser.error("Please provide residue numbers or ranges for each chain specified. Separate selections with a single space.")
     # configure log levels
     config.VERBOSE = True if arguments.verbose else False
     config.QUIET = True if arguments.quiet else False
@@ -268,15 +293,32 @@ def main():
     config.BREAKCOMPOSITE = arguments.breakcomposite
     config.ALTLOC = arguments.altlocation
     config.PEPTIDES = arguments.peptides
+    config.RESIDUES = dict(zip(arguments.peptides, map(residue_list, arguments.residues)))
     config.INTRA = arguments.intra
     config.NOFIX = arguments.nofix
     config.NOFIXFILE = arguments.nofixfile
-    config.NOPDBCANMAP = arguments.nopdbcanmap
+    config.NOPDBCANMAP = bool(arguments.nopdbcanmap or config.INTRA or config.PEPTIDES)
     config.KEEPMOD = arguments.keepmod
     config.DNARECEPTOR = arguments.dnareceptor
     config.OUTPUTFILENAME = arguments.outputfilename
     config.NOHYDRO = arguments.nohydro
     config.MODEL = arguments.model
+
+    try:
+        # add inner quotes for python backend
+        if not arguments.chains:
+            config.CHAINS = None
+        else:
+            import re
+            quoted_input = re.sub(r'(?<!["\'])\b([a-zA-Z0-9_]+)\b(?!["\'])', r'"\1"', arguments.chains)
+            config.CHAINS = ast.literal_eval(quoted_input)
+        print(config.CHAINS)
+        if config.CHAINS and not all(isinstance(c, list) for c in config.CHAINS):
+            raise ValueError("Chains should be specified as a list of lists, e.g., '[[A], [B, C]]'.")
+    except (ValueError, SyntaxError):
+        parser.error("The --chains option must be in the format '[[A], [B, C]]'.")
+
+
     # Make sure we have pymol with --pics and --pymol
     if config.PICS or config.PYMOL:
         try:


=====================================
plip/structure/detection.py
=====================================
@@ -291,7 +291,7 @@ def water_bridges(bs_hba, lig_hba, bs_hbd, lig_hbd, water):
         if not wl.oxy == wd.oxy:
             continue
         # Same water molecule and angle within omega
-        w_angle = vecangle(vector(acc.a.coords, wl.oxy.coords), vector(wl.oxy.coords, don.h.coords))
+        w_angle = vecangle(vector(wl.oxy.coords, acc.a.coords), vector(wl.oxy.coords, don.h.coords))
         if not config.WATER_BRIDGE_OMEGA_MIN < w_angle < config.WATER_BRIDGE_OMEGA_MAX:
             continue
         resnr, reschain, restype = whichresnumber(don.d), whichchain(don.d), whichrestype(don.d)
@@ -309,7 +309,7 @@ def water_bridges(bs_hba, lig_hba, bs_hbd, lig_hbd, water):
         if not wl.oxy == wd.oxy:
             continue
         # Same water molecule and angle within omega
-        w_angle = vecangle(vector(acc.a.coords, wl.oxy.coords), vector(wl.oxy.coords, don.h.coords))
+        w_angle = vecangle(vector(wl.oxy.coords, acc.a.coords), vector(wl.oxy.coords, don.h.coords))
         if not config.WATER_BRIDGE_OMEGA_MIN < w_angle < config.WATER_BRIDGE_OMEGA_MAX:
             continue
         resnr, reschain, restype = whichresnumber(acc.a), whichchain(acc.a), whichrestype(acc.a)


=====================================
plip/structure/preparation.py
=====================================
@@ -14,6 +14,7 @@ from plip.basic.supplemental import cluster_doubles, is_lig, normalize_vector, v
 from plip.basic.supplemental import extract_pdbid, read_pdb, create_folder_if_not_exists, canonicalize
 from plip.basic.supplemental import read, nucleotide_linkage, sort_members_by_importance
 from plip.basic.supplemental import whichchain, whichrestype, whichresnumber, euclidean3d, int32_to_negative
+from plip.basic.supplemental import residue_belongs_to_receptor
 from plip.structure.detection import halogen, pication, water_bridges, metal_complexation
 from plip.structure.detection import hydrophobic_interactions, pistacking, hbonds, saltbridge
 
@@ -243,7 +244,7 @@ class LigandFinder:
         given chain without water
         """
         all_from_chain = [o for o in pybel.ob.OBResidueIter(
-            self.proteincomplex.OBMol) if o.GetChain() == chain]  # All residues from chain
+            self.proteincomplex.OBMol) if o.GetChain() == chain and not self.is_het_residue(o)]  # All residues from chain
         if len(all_from_chain) == 0:
             return None
         else:
@@ -256,7 +257,7 @@ class LigandFinder:
         Returns all non-empty ligands.
         """
 
-        if config.PEPTIDES == [] and config.INTRA is None:
+        if config.PEPTIDES == [] and config.INTRA is None and config.CHAINS is None:
             # Extract small molecule ligands (default)
             ligands = []
 
@@ -284,8 +285,16 @@ class LigandFinder:
         else:
             # Extract peptides from given chains
             self.water = [o for o in pybel.ob.OBResidueIter(self.proteincomplex.OBMol) if o.GetResidueProperty(9)]
-            if config.PEPTIDES:
+            if config.PEPTIDES and not config.CHAINS:
                 peptide_ligands = [self.getpeptides(chain) for chain in config.PEPTIDES]
+
+            #Todo: Validate change here... Do we want to combine multiple chains to a single ligand?
+            # if yes can be easily added to the getpeptides function by flatten the resulting list - Philipp
+            elif config.CHAINS:
+                # chains is defined as list of list e.g. [['A'], ['B', 'C']] in which second list contains the
+                # ligand chains and the first one should be the receptor
+                peptide_ligands = [self.getpeptides(chain) for chain in config.CHAINS[1]]
+
             elif config.INTRA is not None:
                 peptide_ligands = [self.getpeptides(config.INTRA), ]
 
@@ -304,7 +313,7 @@ class LigandFinder:
         names = [x[0] for x in members]
         longname = '-'.join([x[0] for x in members])
 
-        if config.PEPTIDES:
+        if config.PEPTIDES or config.CHAINS:
             ligtype = 'PEPTIDE'
         elif config.INTRA is not None:
             ligtype = 'INTRA'
@@ -328,26 +337,20 @@ class LigandFinder:
 
         # create a new ligand molecule representing any k-mer linked structures
         lig = pybel.ob.OBMol()
-        neighbours = dict()
-        for obatom in hetatoms.values():  # iterate over atom objects
-            idx = obatom.GetIdx()
-            lig.AddAtom(obatom)
-            # ids of all neighbours of obatom
-            neighbours[idx] = set([neighbour_atom.GetIdx() for neighbour_atom
-                                   in pybel.ob.OBAtomAtomIter(obatom)]) & set(hetatoms.keys())
-        logger.debug(f'atom neighbours mapped')
+        # create bit vector with 1 for all ligand atoms and 0 for all others
+        atomsBitVec = pybel.ob.OBBitVec(self.proteincomplex.OBMol.NumAtoms())
+        for atomidx in hetatoms.keys():
+            atomsBitVec.SetBitOn(atomidx)
+        # safe substructure defined by those atoms to previously empty lig
+        self.proteincomplex.OBMol.CopySubstructure(lig, atomsBitVec, None, 0)
 
         ##############################################################
         # map the old atom idx of OBMol to the new idx of the ligand #
         ##############################################################
 
-        newidx = dict(zip(hetatoms.keys(), [obatom.GetIdx() for obatom in pybel.ob.OBMolAtomIter(lig)]))
+        newidx = dict(zip(sorted(hetatoms.keys()), [obatom.GetIdx() for obatom in pybel.ob.OBMolAtomIter(lig)]))
         mapold = dict(zip(newidx.values(), newidx))
-        # copy the bonds
-        for obatom in hetatoms:
-            for neighbour_atom in neighbours[obatom]:
-                bond = hetatoms[obatom].GetBond(hetatoms[neighbour_atom])
-                lig.AddBond(newidx[obatom], newidx[neighbour_atom], bond.GetBondOrder())
+
         lig = pybel.Molecule(lig)
 
         # For kmers, the representative ids are chosen (first residue of kmer)
@@ -515,11 +518,12 @@ class Mol:
         """Find all possible hydrogen bond acceptors"""
         data = namedtuple('hbondacceptor', 'a a_orig_atom a_orig_idx type')
         a_set = []
-        for atom in filter(lambda at: at.OBAtom.IsHbondAcceptor(), all_atoms):
+        for atom in all_atoms:
             if atom.atomicnum not in [9, 17, 35, 53] and atom.idx not in self.altconf:  # Exclude halogen atoms
                 a_orig_idx = self.Mapper.mapid(atom.idx, mtype=self.mtype, bsid=self.bsid)
                 a_orig_atom = self.Mapper.id_to_atom(a_orig_idx)
-                a_set.append(data(a=atom, a_orig_atom=a_orig_atom, a_orig_idx=a_orig_idx, type='regular'))
+                if a_orig_atom.OBAtom.IsHbondAcceptor():
+                    a_set.append(data(a=atom, a_orig_atom=a_orig_atom, a_orig_idx=a_orig_idx, type='regular'))
         a_set = sorted(a_set, key=lambda x: x.a_orig_idx)
         return a_set
 
@@ -528,13 +532,11 @@ class Mol:
         donor_pairs = []
         data = namedtuple('hbonddonor', 'd d_orig_atom d_orig_idx h type')
         for donor in [a for a in all_atoms if a.OBAtom.IsHbondDonor() and a.idx not in self.altconf]:
-            in_ring = False
-            if not in_ring:
-                for adj_atom in [a for a in pybel.ob.OBAtomAtomIter(donor.OBAtom) if a.IsHbondDonorH()]:
-                    d_orig_idx = self.Mapper.mapid(donor.idx, mtype=self.mtype, bsid=self.bsid)
-                    d_orig_atom = self.Mapper.id_to_atom(d_orig_idx)
-                    donor_pairs.append(data(d=donor, d_orig_atom=d_orig_atom, d_orig_idx=d_orig_idx,
-                                            h=pybel.Atom(adj_atom), type='regular'))
+            for adj_atom in [a for a in pybel.ob.OBAtomAtomIter(donor.OBAtom) if a.IsHbondDonorH()]:
+                d_orig_idx = self.Mapper.mapid(donor.idx, mtype=self.mtype, bsid=self.bsid)
+                d_orig_atom = self.Mapper.id_to_atom(d_orig_idx)
+                donor_pairs.append(data(d=donor, d_orig_atom=d_orig_atom, d_orig_idx=d_orig_idx,
+                                        h=pybel.Atom(adj_atom), type='regular'))
         for carbon in hydroph_atoms:
             for adj_atom in [a for a in pybel.ob.OBAtomAtomIter(carbon.atom.OBAtom) if a.GetAtomicNum() == 1]:
                 d_orig_idx = self.Mapper.mapid(carbon.atom.idx, mtype=self.mtype, bsid=self.bsid)
@@ -631,10 +633,11 @@ class PLInteraction:
 
         self.pistacking = pistacking(self.bindingsite.rings, self.ligand.rings)
 
-        self.all_pi_cation_laro = pication(self.ligand.rings, self.bindingsite.get_pos_charged(), True)
-        self.pication_paro = pication(self.bindingsite.rings, self.ligand.get_pos_charged(), False)
+        self.all_pication_laro = pication(self.ligand.rings, self.bindingsite.get_pos_charged(), True)
+        self.all_pication_paro = pication(self.bindingsite.rings, self.ligand.get_pos_charged(), False)
 
-        self.pication_laro = self.refine_pi_cation_laro(self.all_pi_cation_laro, self.pistacking)
+        self.pication_laro = self.refine_pication(self.all_pication_laro, self.pistacking)
+        self.pication_paro = self.refine_pication(self.all_pication_paro, self.pistacking)
 
         self.all_hydrophobic_contacts = hydrophobic_interactions(self.bindingsite.get_hydrophobic_atoms(),
                                                                  self.ligand.get_hydrophobic_atoms())
@@ -748,47 +751,59 @@ class PLInteraction:
                     sel2[(h.ligatom.idx, h.resnr)] = h
         hydroph = [h for h in sel2.values()]
         hydroph_final = []
-        bsclust = {}
         #  3. If a protein atom interacts with several neighboring ligand atoms, just keep the one with the closest dist
-        for h in hydroph:
-            if h.bsatom.idx not in bsclust:
-                bsclust[h.bsatom.idx] = [h, ]
-            else:
-                bsclust[h.bsatom.idx].append(h)
-
-        idx_to_h = {}
-        for bs in [a for a in bsclust if len(bsclust[a]) == 1]:
-            hydroph_final.append(bsclust[bs][0])
-
-        # A list of tuples with the idx of an atom and one of its neighbours is created
-        for bs in [a for a in bsclust if not len(bsclust[a]) == 1]:
-            tuples = []
-            all_idx = [i.ligatom.idx for i in bsclust[bs]]
-            for b in bsclust[bs]:
-                idx = b.ligatom.idx
-                neigh = [na for na in pybel.ob.OBAtomAtomIter(b.ligatom.OBAtom)]
-                for n in neigh:
-                    n_idx = n.GetIdx()
-                    if n_idx in all_idx:
-                        if n_idx < idx:
-                            tuples.append((n_idx, idx))
-                        else:
-                            tuples.append((idx, n_idx))
-                        idx_to_h[idx] = b
-
-            tuples = list(set(tuples))
-            tuples = sorted(tuples, key=itemgetter(1))
-            clusters = cluster_doubles(tuples)  # Cluster connected atoms (i.e. find hydrophobic patches)
-
-            for cluster in clusters:
-                min_dist = float('inf')
-                min_h = None
-                for atm_idx in cluster:
-                    h = idx_to_h[atm_idx]
-                    if h.distance < min_dist:
-                        min_dist = h.distance
-                        min_h = h
-                hydroph_final.append(min_h)
+        if config.PEPTIDES or config.INTRA or config.CHAINS:
+            # the ligand also consists of amino acid residues, repeat step 2 just the other way around
+            sel3 = {}
+            for h in hydroph:
+                if not (h.bsatom.idx, h.resnr_l) in sel3:
+                    sel3[(h.bsatom.idx, h.resnr_l)] = h
+                else:
+                    if sel3[(h.bsatom.idx, h.resnr_l)].distance > h.distance:
+                        sel3[(h.bsatom.idx, h.resnr_l)] = h
+            hydroph_final = [h for h in sel3.values()]
+        else:
+            # for other ligands find hydrophobic patches for which only the shortest interactions is reported
+            bsclust = {}
+            for h in hydroph:
+                if h.bsatom.idx not in bsclust:
+                    bsclust[h.bsatom.idx] = [h, ]
+                else:
+                    bsclust[h.bsatom.idx].append(h)
+
+            idx_to_h = {}
+            for bs in [a for a in bsclust if len(bsclust[a]) == 1]:
+                hydroph_final.append(bsclust[bs][0])
+
+            # A list of tuples with the idx of an atom and one of its neighbours is created
+            for bs in [a for a in bsclust if not len(bsclust[a]) == 1]:
+                tuples = []
+                all_idx = [i.ligatom.idx for i in bsclust[bs]]
+                for b in bsclust[bs]:
+                    idx = b.ligatom.idx
+                    neigh = [na for na in pybel.ob.OBAtomAtomIter(b.ligatom.OBAtom)]
+                    for n in neigh:
+                        n_idx = n.GetIdx()
+                        if n_idx in all_idx:
+                            if n_idx < idx:
+                                tuples.append((n_idx, idx))
+                            else:
+                                tuples.append((idx, n_idx))
+                            idx_to_h[idx] = b
+
+                tuples = list(set(tuples))
+                tuples = sorted(tuples, key=itemgetter(1))
+                clusters = cluster_doubles(tuples)  # Cluster connected atoms (i.e. find hydrophobic patches)
+
+                for cluster in clusters:
+                    min_dist = float('inf')
+                    min_h = None
+                    for atm_idx in cluster:
+                        h = idx_to_h[atm_idx]
+                        if h.distance < min_dist:
+                            min_dist = h.distance
+                            min_h = h
+                    hydroph_final.append(min_h)
         before, reduced = len(all_h), len(hydroph_final)
         if not before == 0 and not before == reduced:
             logger.info(f'reduced number of hydrophobic contacts from {before} to {reduced}')
@@ -849,7 +864,7 @@ class PLInteraction:
         return [hb[1] for hb in second_set.values()]
 
     @staticmethod
-    def refine_pi_cation_laro(all_picat, stacks):
+    def refine_pication(all_picat, stacks):
         """Just important for constellations with histidine involved. If the histidine ring is positioned in stacking
         position to an aromatic ring in the ligand, there is in most cases stacking and pi-cation interaction reported
         as histidine also carries a positive charge in the ring. For such cases, only report stacking.
@@ -857,9 +872,15 @@ class PLInteraction:
         i_set = []
         for picat in all_picat:
             exclude = False
-            for stack in stacks:
-                if whichrestype(stack.proteinring.atoms[0]) == 'HIS' and picat.ring.obj == stack.ligandring.obj:
-                    exclude = True
+            if picat.restype == 'HIS':
+                for stack in stacks:
+                    if whichresnumber(stack.proteinring.atoms[0]) == picat.resnr and picat.ring.obj == stack.ligandring.obj:
+                        exclude = True
+            # HIS could also be on ligand side for protein-protein interactions
+            if picat.restype_l == 'HIS':
+                for stack in stacks:
+                    if whichresnumber(stack.ligandring.atoms[0]) == picat.resnr_l and picat.ring.obj == stack.proteinring.obj:
+                        exclude = True
             if not exclude:
                 i_set.append(picat)
         return i_set
@@ -868,18 +889,24 @@ class PLInteraction:
     def refine_water_bridges(wbridges, hbonds_ldon, hbonds_pdon):
         """A donor atom already forming a hydrogen bond is not allowed to form a water bridge. Each water molecule
         can only be donor for two water bridges, selecting the constellation with the omega angle closest to 110 deg."""
-        donor_atoms_hbonds = [hb.d.idx for hb in hbonds_ldon + hbonds_pdon]
+        donor_atoms_hbonds = [hb.d_orig_idx for hb in hbonds_ldon + hbonds_pdon]
         wb_dict = {}
         wb_dict2 = {}
         omega = 110.0
 
         # Just one hydrogen bond per donor atom
-        for wbridge in [wb for wb in wbridges if wb.d.idx not in donor_atoms_hbonds]:
-            if (wbridge.water.idx, wbridge.a.idx) not in wb_dict:
-                wb_dict[(wbridge.water.idx, wbridge.a.idx)] = wbridge
+        # (TODO: looks wrong, what about donor atoms with multiple H? Does PLIP cover that somehow?
+        # and only one donor per water-acceptor combination
+        # (each water-acceptor combo only kept once after this loop,
+        # the one with angle at water closest to omega)
+        for wbridge in [wb for wb in wbridges if wb.d_orig_idx not in donor_atoms_hbonds]:
+            if (wbridge.water_orig_idx, wbridge.a_orig_idx) not in wb_dict:
+                wb_dict[(wbridge.water_orig_idx, wbridge.a_orig_idx)] = wbridge
             else:
-                if abs(omega - wb_dict[(wbridge.water.idx, wbridge.a.idx)].w_angle) < abs(omega - wbridge.w_angle):
-                    wb_dict[(wbridge.water.idx, wbridge.a.idx)] = wbridge
+                if abs(omega - wb_dict[(wbridge.water_orig_idx, wbridge.a_orig_idx)].w_angle) > abs(omega - wbridge.w_angle):
+                    wb_dict[(wbridge.water_orig_idx, wbridge.a_orig_idx)] = wbridge
+        # Just two hydrogen bonds per water molecule
+        # only the two water-acceptor combos with angle closest to omega are kept
         for wb_tuple in wb_dict:
             water, acceptor = wb_tuple
             if water not in wb_dict2:
@@ -888,7 +915,7 @@ class PLInteraction:
                 wb_dict2[water].append((abs(omega - wb_dict[wb_tuple].w_angle), wb_dict[wb_tuple]))
                 wb_dict2[water] = sorted(wb_dict2[water], key=lambda x: x[0])
             else:
-                if wb_dict2[water][1][0] < abs(omega - wb_dict[wb_tuple].w_angle):
+                if wb_dict2[water][1][0] > abs(omega - wb_dict[wb_tuple].w_angle):
                     wb_dict2[water] = [wb_dict2[water][0], (wb_dict[wb_tuple].w_angle, wb_dict[wb_tuple])]
 
         filtered_wb = []
@@ -933,7 +960,7 @@ class BindingSite(Mol):
         data = namedtuple('pcharge', 'atoms atoms_orig_idx type center restype resnr reschain')
         a_set = []
         # Iterate through all residue, exclude those in chains defined as peptides
-        for res in [r for r in pybel.ob.OBResidueIter(mol.OBMol) if not r.GetChain() in config.PEPTIDES]:
+        for res in [r for r in pybel.ob.OBResidueIter(mol.OBMol) if residue_belongs_to_receptor(r, config)]:
             if config.INTRA is not None:
                 if res.GetChain() != config.INTRA:
                     continue
@@ -974,7 +1001,7 @@ class BindingSite(Mol):
                         a_contributing.append(pybel.Atom(a))
                         a_contributing_orig_idx.append(self.Mapper.mapid(a.GetIdx(), mtype='protein'))
                 if not len(a_contributing) == 0:
-                    a_set.append(data(atoms=a_contributing,atoms_orig_idx=a_contributing_orig_idx, type='negative', 
+                    a_set.append(data(atoms=a_contributing,atoms_orig_idx=a_contributing_orig_idx, type='negative',
                                       center=centroid([ac.coords for ac in a_contributing]), restype=res.GetName(),
                                       resnr=res.GetNum(),
                                       reschain=res.GetChain()))
@@ -1032,7 +1059,8 @@ class Ligand(Mol):
         self.type = ligand.type
         self.complex = cclass
         self.molecule = ligand.mol  # Pybel Molecule
-        self.smiles = self.molecule.write(format='can')  # SMILES String
+        # get canonical SMILES String, but not for peptide ligand (tend to be too long -> openBabel crashes)
+        self.smiles = "" if (config.INTRA or config.PEPTIDES or config.CHAINS) else self.molecule.write(format='can')
         self.inchikey = self.molecule.write(format='inchikey')
         self.can_to_pdb = ligand.can_to_pdb
         if not len(self.smiles) == 0:
@@ -1173,74 +1201,110 @@ class Ligand(Mol):
         """
         data = namedtuple('lcharge', 'atoms orig_atoms atoms_orig_idx type center fgroup')
         a_set = []
-        for a in all_atoms:
-            a_orig_idx = self.Mapper.mapid(a.idx, mtype=self.mtype, bsid=self.bsid)
-            a_orig = self.Mapper.id_to_atom(a_orig_idx)
-            if self.is_functional_group(a, 'quartamine'):
-                a_set.append(data(atoms=[a, ], orig_atoms=[a_orig, ], atoms_orig_idx=[a_orig_idx, ], type='positive',
-                                  center=list(a.coords), fgroup='quartamine'))
-            elif self.is_functional_group(a, 'tertamine'):
-                a_set.append(data(atoms=[a, ], orig_atoms=[a_orig, ], atoms_orig_idx=[a_orig_idx, ], type='positive',
-                                  center=list(a.coords),
-                                  fgroup='tertamine'))
-            if self.is_functional_group(a, 'sulfonium'):
-                a_set.append(data(atoms=[a, ], orig_atoms=[a_orig, ], atoms_orig_idx=[a_orig_idx, ], type='positive',
-                                  center=list(a.coords),
-                                  fgroup='sulfonium'))
-            if self.is_functional_group(a, 'phosphate'):
-                a_contributing = [a, ]
-                a_contributing_orig_idx = [a_orig_idx, ]
-                [a_contributing.append(pybel.Atom(neighbor)) for neighbor in pybel.ob.OBAtomAtomIter(a.OBAtom)]
-                [a_contributing_orig_idx.append(self.Mapper.mapid(neighbor.idx, mtype=self.mtype, bsid=self.bsid))
-                 for neighbor in a_contributing]
-                orig_contributing = [self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx]
-                a_set.append(
-                    data(atoms=a_contributing, orig_atoms=orig_contributing, atoms_orig_idx=a_contributing_orig_idx,
-                         type='negative',
-                         center=a.coords, fgroup='phosphate'))
-            if self.is_functional_group(a, 'sulfonicacid'):
-                a_contributing = [a, ]
-                a_contributing_orig_idx = [a_orig_idx, ]
-                [a_contributing.append(pybel.Atom(neighbor)) for neighbor in pybel.ob.OBAtomAtomIter(a.OBAtom) if
-                 neighbor.GetAtomicNum() == 8]
-                [a_contributing_orig_idx.append(self.Mapper.mapid(neighbor.idx, mtype=self.mtype, bsid=self.bsid))
-                 for neighbor in a_contributing]
-                orig_contributing = [self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx]
-                a_set.append(
-                    data(atoms=a_contributing, orig_atoms=orig_contributing, atoms_orig_idx=a_contributing_orig_idx,
-                         type='negative',
-                         center=a.coords, fgroup='sulfonicacid'))
-            elif self.is_functional_group(a, 'sulfate'):
-                a_contributing = [a, ]
-                a_contributing_orig_idx = [a_orig_idx, ]
-                [a_contributing_orig_idx.append(self.Mapper.mapid(neighbor.idx, mtype=self.mtype, bsid=self.bsid))
-                 for neighbor in a_contributing]
-                [a_contributing.append(pybel.Atom(neighbor)) for neighbor in pybel.ob.OBAtomAtomIter(a.OBAtom)]
-                orig_contributing = [self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx]
-                a_set.append(
-                    data(atoms=a_contributing, orig_atoms=orig_contributing, atoms_orig_idx=a_contributing_orig_idx,
-                         type='negative',
-                         center=a.coords, fgroup='sulfate'))
-            if self.is_functional_group(a, 'carboxylate'):
-                a_contributing = [pybel.Atom(neighbor) for neighbor in pybel.ob.OBAtomAtomIter(a.OBAtom)
-                                  if neighbor.GetAtomicNum() == 8]
-                a_contributing_orig_idx = [self.Mapper.mapid(neighbor.idx, mtype=self.mtype, bsid=self.bsid)
-                                           for neighbor in a_contributing]
-                orig_contributing = [self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx]
-                a_set.append(
-                    data(atoms=a_contributing, orig_atoms=orig_contributing, atoms_orig_idx=a_contributing_orig_idx,
-                         type='negative',
-                         center=centroid([a.coords for a in a_contributing]), fgroup='carboxylate'))
-            elif self.is_functional_group(a, 'guanidine'):
-                a_contributing = [pybel.Atom(neighbor) for neighbor in pybel.ob.OBAtomAtomIter(a.OBAtom)
-                                  if neighbor.GetAtomicNum() == 7]
-                a_contributing_orig_idx = [self.Mapper.mapid(neighbor.idx, mtype=self.mtype, bsid=self.bsid)
-                                           for neighbor in a_contributing]
-                orig_contributing = [self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx]
-                a_set.append(
-                    data(atoms=a_contributing, orig_atoms=orig_contributing, atoms_orig_idx=a_contributing_orig_idx,
-                         type='positive',
-                         center=a.coords, fgroup='guanidine'))
+        if not (config.INTRA or config.PEPTIDES or config.CHAINS):
+            for a in all_atoms:
+                a_orig_idx = self.Mapper.mapid(a.idx, mtype=self.mtype, bsid=self.bsid)
+                a_orig = self.Mapper.id_to_atom(a_orig_idx)
+                if self.is_functional_group(a, 'quartamine'):
+                    a_set.append(data(atoms=[a, ], orig_atoms=[a_orig, ], atoms_orig_idx=[a_orig_idx, ], type='positive',
+                                      center=list(a.coords), fgroup='quartamine'))
+                elif self.is_functional_group(a, 'tertamine'):
+                    a_set.append(data(atoms=[a, ], orig_atoms=[a_orig, ], atoms_orig_idx=[a_orig_idx, ], type='positive',
+                                      center=list(a.coords),
+                                      fgroup='tertamine'))
+                if self.is_functional_group(a, 'sulfonium'):
+                    a_set.append(data(atoms=[a, ], orig_atoms=[a_orig, ], atoms_orig_idx=[a_orig_idx, ], type='positive',
+                                      center=list(a.coords),
+                                      fgroup='sulfonium'))
+                if self.is_functional_group(a, 'phosphate'):
+                    a_contributing = [a, ]
+                    a_contributing_orig_idx = [a_orig_idx, ]
+                    [a_contributing.append(pybel.Atom(neighbor)) for neighbor in pybel.ob.OBAtomAtomIter(a.OBAtom)]
+                    [a_contributing_orig_idx.append(self.Mapper.mapid(neighbor.idx, mtype=self.mtype, bsid=self.bsid))
+                     for neighbor in a_contributing]
+                    orig_contributing = [self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx]
+                    a_set.append(
+                        data(atoms=a_contributing, orig_atoms=orig_contributing, atoms_orig_idx=a_contributing_orig_idx,
+                             type='negative',
+                             center=a.coords, fgroup='phosphate'))
+                if self.is_functional_group(a, 'sulfonicacid'):
+                    a_contributing = [a, ]
+                    a_contributing_orig_idx = [a_orig_idx, ]
+                    [a_contributing.append(pybel.Atom(neighbor)) for neighbor in pybel.ob.OBAtomAtomIter(a.OBAtom) if
+                     neighbor.GetAtomicNum() == 8]
+                    [a_contributing_orig_idx.append(self.Mapper.mapid(neighbor.idx, mtype=self.mtype, bsid=self.bsid))
+                     for neighbor in a_contributing]
+                    orig_contributing = [self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx]
+                    a_set.append(
+                        data(atoms=a_contributing, orig_atoms=orig_contributing, atoms_orig_idx=a_contributing_orig_idx,
+                             type='negative',
+                             center=a.coords, fgroup='sulfonicacid'))
+                elif self.is_functional_group(a, 'sulfate'):
+                    a_contributing = [a, ]
+                    a_contributing_orig_idx = [a_orig_idx, ]
+                    [a_contributing_orig_idx.append(self.Mapper.mapid(neighbor.idx, mtype=self.mtype, bsid=self.bsid))
+                     for neighbor in a_contributing]
+                    [a_contributing.append(pybel.Atom(neighbor)) for neighbor in pybel.ob.OBAtomAtomIter(a.OBAtom)]
+                    orig_contributing = [self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx]
+                    a_set.append(
+                        data(atoms=a_contributing, orig_atoms=orig_contributing, atoms_orig_idx=a_contributing_orig_idx,
+                             type='negative',
+                             center=a.coords, fgroup='sulfate'))
+                if self.is_functional_group(a, 'carboxylate'):
+                    a_contributing = [pybel.Atom(neighbor) for neighbor in pybel.ob.OBAtomAtomIter(a.OBAtom)
+                                      if neighbor.GetAtomicNum() == 8]
+                    a_contributing_orig_idx = [self.Mapper.mapid(neighbor.idx, mtype=self.mtype, bsid=self.bsid)
+                                               for neighbor in a_contributing]
+                    orig_contributing = [self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx]
+                    a_set.append(
+                        data(atoms=a_contributing, orig_atoms=orig_contributing, atoms_orig_idx=a_contributing_orig_idx,
+                             type='negative',
+                             center=centroid([a.coords for a in a_contributing]), fgroup='carboxylate'))
+                elif self.is_functional_group(a, 'guanidine'):
+                    a_contributing = [pybel.Atom(neighbor) for neighbor in pybel.ob.OBAtomAtomIter(a.OBAtom)
+                                      if neighbor.GetAtomicNum() == 7]
+                    a_contributing_orig_idx = [self.Mapper.mapid(neighbor.idx, mtype=self.mtype, bsid=self.bsid)
+                                               for neighbor in a_contributing]
+                    orig_contributing = [self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx]
+                    a_set.append(
+                        data(atoms=a_contributing, orig_atoms=orig_contributing, atoms_orig_idx=a_contributing_orig_idx,
+                             type='positive',
+                             center=a.coords, fgroup='guanidine'))
+        else:
+            """We have peptide/protein chain as ligand"""
+            """Looks for positive charges in arginine, histidine or lysine, for negative in aspartic and glutamic acid."""
+            for res in pybel.ob.OBResidueIter(self.molecule.OBMol):
+                if config.INTRA is not None:
+                    if res.GetChain() != config.INTRA:
+                        continue
+                a_contributing = []
+                a_contributing_orig_idx = []
+                if res.GetName() in ('ARG', 'HIS', 'LYS'):  # Arginine, Histidine or Lysine have charged sidechains
+                    for a in pybel.ob.OBResidueAtomIter(res):
+                        if a.GetType().startswith('N') and res.GetAtomProperty(a, 8) \
+                                and not self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid) in self.altconf:
+                            a_contributing.append(pybel.Atom(a))
+                            a_contributing_orig_idx.append(self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid))
+                    if not len(a_contributing) == 0:
+                        a_set.append(data(atoms=a_contributing,
+                                          orig_atoms=[self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx],
+                                          atoms_orig_idx=a_contributing_orig_idx,
+                                          type='positive',
+                                          center=centroid([ac.coords for ac in a_contributing]),
+                                          fgroup=res.GetName()+str(res.GetNum())+res.GetChain()))
+                if res.GetName() in ('GLU', 'ASP'):  # Aspartic or Glutamic Acid
+                    for a in pybel.ob.OBResidueAtomIter(res):
+                        if a.GetType().startswith('O') and res.GetAtomProperty(a, 8) \
+                                and not self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid) in self.altconf:
+                            a_contributing.append(pybel.Atom(a))
+                            a_contributing_orig_idx.append(self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid))
+                    if not len(a_contributing) == 0:
+                        a_set.append(data(atoms=a_contributing,
+                                          orig_atoms=[self.Mapper.id_to_atom(idx) for idx in a_contributing_orig_idx],
+                                          atoms_orig_idx=a_contributing_orig_idx,
+                                          type='negative',
+                                          center=centroid([ac.coords for ac in a_contributing]),
+                                          fgroup=res.GetName()+str(res.GetNum())+res.GetChain()))
         return a_set
 
     def find_metal_binding(self, lig_atoms, water_oxygens):
@@ -1512,9 +1576,9 @@ class PDBComplex:
         rescentroid = centroid([(atm.x(), atm.y(), atm.z()) for atm in pybel.ob.OBResidueAtomIter(res)])
         # Check geometry
         near_enough = True if euclidean3d(rescentroid, ligcentroid) < cutoff else False
-        # Check chain membership
-        restricted_chain = True if res.GetChain() in config.PEPTIDES else False
-        return (near_enough and not restricted_chain)
+        #Todo: Test if properly working
+        # Add restriction via chains flag
+        return near_enough and residue_belongs_to_receptor(res, config)
 
     def get_atom(self, idx):
         return self.atoms[idx]
@@ -1525,4 +1589,4 @@ class PDBComplex:
 
     @output_path.setter
     def output_path(self, path):
-        self._output_path = tilde_expansion(path)
+        self._output_path = tilde_expansion(path)
\ No newline at end of file


=====================================
plip/test/test_literature_validated.py
=====================================
@@ -232,6 +232,9 @@ class LiteratureValidatedTest(unittest.TestCase):
         # Water bridges to Lys307, Arg309 and 111 from phosphate groups
         waterbridges = {wb.resnr for wb in s.water_bridges}
         # Water bridge to 60 not found due to prioritization
+        # Philipp: No interaction to 111 due to prioritization/omega angle compute correction
+        # --> Now two interactions with ASN92 which were before only one that is not checked here
+        # self.assertTrue({307, 309}.issubset(waterbridges))
         self.assertTrue({307, 309, 111}.issubset(waterbridges))
         # pi-stacking interaction with Phe209
         pistackres = {pistack.resnr for pistack in s.pistacking}
@@ -312,6 +315,8 @@ class LiteratureValidatedTest(unittest.TestCase):
                 tmpmol.characterize_complex(ligand)
         s = tmpmol.interaction_sets[bsid]
         # Hydrogen bonds to Ala609
+        # Philipp: In this case the protein is the acceptor
+        # hbonds = {hbond.resnr for hbond in s.hbonds_pdon + s.hbonds_ldon}
         hbonds = {hbond.resnr for hbond in s.hbonds_pdon}
         self.assertTrue({609}.issubset(hbonds))
         # Saltbridge to Asp513
@@ -515,12 +520,16 @@ class LiteratureValidatedTest(unittest.TestCase):
         s = tmpmol.interaction_sets[bsid]
         #@todo Publication show hydrogen bond interactions for Gly219
         # Hydrogen bonds to Ser190, Ser195, Gly219 and Asp189
+
+        #Philipp: Only Hbont to 189 + 195 are possible. Else, donors would be included in multiple interactions
         hbonds = {hbond.resnr for hbond in s.hbonds_pdon+s.hbonds_ldon}
-        self.assertTrue({189, 190, 195}.issubset(hbonds))
+        #self.assertTrue({189, 190, 195}.issubset(hbonds))
+        self.assertTrue({189, 195}.issubset(hbonds))
         # Water bridges to Ser190 and Val227
         # Water bridge to 190 not detected due to prioritization
+        # Philipp: Listen to the comment above
         waterbridges = {wb.resnr for wb in s.water_bridges}
-        self.assertTrue({227}.issubset(waterbridges))
+        #self.assertTrue({227}.issubset(waterbridges))
         # hydrophobic interaction of Leu99
         hydrophobics = {hydrophobic.resnr for hydrophobic in s.all_hydrophobic_contacts}
         self.assertTrue({99}.issubset(hydrophobics))
@@ -587,6 +596,8 @@ class LiteratureValidatedTest(unittest.TestCase):
         # Water bridges to Trp137
         # res nr 52 mentioned in alternative conformation not considered
         waterbridges = {wb.resnr for wb in s.water_bridges}
+        # Philipp: Now detects two wbridges with ARG57 because the omega angles are better
+        #self.assertTrue({57}.issubset(waterbridges))
         self.assertTrue({322}.issubset(waterbridges))
         # pi-stacking interaction with Trp384, Trp137 and Trp52
         pistackres = {pistack.resnr for pistack in s.pistacking}
@@ -701,6 +712,8 @@ class LiteratureValidatedTest(unittest.TestCase):
         # Water bridges
         waterbridges = {str(wb.resnr)+wb.reschain for wb in s.water_bridges}
         # #@todo Water bridge with 50B not detected
+        # Philipp: Water bridge with 50A not detected -> Angle prio leads to two wbridges with 50B
+        # self.assertTrue({'50B'}.issubset(waterbridges))
         self.assertTrue({'50A'}.issubset(waterbridges))  # Bridging Ile-B50 and Ile-A50 with ligand
         # pi-cation Interactions
         picat = {pication.resnr for pication in s.pication_laro}
@@ -752,4 +765,6 @@ class LiteratureValidatedTest(unittest.TestCase):
         # Water bridges
         waterbridges = {str(wb.resnr)+wb.reschain for wb in s.water_bridges}
         # Waterbridge with Gly27 is detected instead of Ala28/Asp29
+        # Philipp: Waterbridge with 50B not detected rather two towards 50A
+        # self.assertTrue({'50B', '29A'}.issubset(waterbridges))
         self.assertTrue({'50A', '50B', '29A'}.issubset(waterbridges))


=====================================
plip/visualization/visualize.py
=====================================
@@ -19,7 +19,7 @@ def visualize_in_pymol(plcomplex):
     pdbid = plcomplex.pdbid
     lig_members = plcomplex.lig_members
     chain = plcomplex.chain
-    if config.PEPTIDES:
+    if config.PEPTIDES or config.CHAINS:
         vis.ligname = 'PeptideChain%s' % plcomplex.chain
     if config.INTRA is not None:
         vis.ligname = 'Intra%s' % plcomplex.chain
@@ -45,7 +45,10 @@ def visualize_in_pymol(plcomplex):
     cmd.set_name(current_name, pdbid)
     cmd.hide('everything', 'all')
     if config.PEPTIDES:
-        cmd.select(ligname, 'chain %s and not resn HOH' % plcomplex.chain)
+        if plcomplex.chain in config.RESIDUES.keys():
+            cmd.select(ligname, 'chain %s and not resn HOH and resi %s' % (plcomplex.chain, "+".join(map(str, config.RESIDUES[plcomplex.chain]))))
+        else:
+            cmd.select(ligname, 'chain %s and not resn HOH' % plcomplex.chain)
     else:
         cmd.select(ligname, 'resn %s and chain %s and resi %s*' % (hetid, chain, plcomplex.position))
     logger.debug(f'selecting ligand for PDBID {pdbid} and ligand name {ligname}')
@@ -97,7 +100,7 @@ def visualize_in_pymol(plcomplex):
         cmd.hide('cartoon', '%sLines' % plcomplex.pdbid)
         cmd.show('lines', '%sLines' % plcomplex.pdbid)
 
-    if config.PEPTIDES:
+    if config.PEPTIDES or config.CHAINS:
         filename = "%s_PeptideChain%s" % (pdbid.upper(), plcomplex.chain)
         if config.PYMOL:
             vis.save_session(config.OUTPATH, override=filename)


=====================================
setup.py
=====================================
@@ -1,32 +1,103 @@
 from setuptools import setup, find_packages
+from setuptools.command.build_ext import build_ext
+from setuptools.command.install import install
+from distutils.command.build import build
 
 from plip.basic import config
 
+def install_pkg_via_pip(package):
+    import sys
+    import subprocess
+    subprocess.check_call([sys.executable, '-m', 'pip', 'install', package])
+
+def build_ob_py_bindings():
+    """ Fix openbabel issue (https://github.com/openbabel/openbabel/issues/2408) manually
+      - Download openbabel bindings from pypi
+      - extract to a tmpdir
+      - fix version in $tmpdir/openbabel-3.1.1.1/openbabel/__init__.py to have 2 dot version only
+      - install the fixed version with setuptools/pip
+      - cleanup the tmpdir
+    """
+    import requests
+    import tarfile
+    import tempfile
+    import shutil
+    import fileinput
+    import subprocess
+    import sys
+
+    openbabel_pypi_url='https://files.pythonhosted.org/packages/9d/3f/f08f5d1422d74ed0e1e612870b343bfcc26313bdf9efec9165c3ea4b3ae2/openbabel-3.1.1.1.tar.gz'
+
+    print (f"Downloading openbabel package from : {openbabel_pypi_url}")
+    obtar=requests.get(openbabel_pypi_url)
+    obtmpdir = tempfile.mkdtemp()
+    obtmp = obtmpdir+'/openbabel-3.1.1.1.tar'
+    open(obtmp,'wb').write(obtar.content)
+    print(f"Saving openbabel tar.gz to {obtmpdir}")
+    versfile = obtmpdir+'/openbabel-3.1.1.1/openbabel/__init__.py'
+
+    with tarfile.open(obtmp,mode='r') as tf:
+        tf.extractall(obtmpdir)
+
+    print ('Fix versions: s/3.1.1.1/3.1.1/ to make StrictVersion() in openbabel\'s setup.py happy')
+    print ('See https://github.com/openbabel/openbabel/issues/2408 for more details')
+    with fileinput.input(files=versfile,inplace=True) as f:
+        for line in f:
+            op = line.replace('__version__ = "3.1.1.1"', '__version__ = "3.1.1"')
+            print(op, end='')
+
+    install_pkg_via_pip(obtmpdir+'/openbabel-3.1.1.1/')
+    print (f"Cleanup tmpdir: {obtmpdir}")
+    shutil.rmtree(obtmpdir)
+
+class CustomBuild(build):
+    """Ensure build_ext runs first in build command."""
+    def run(self):
+        self.run_command('build_ext')
+        build.run(self)
+
+class CustomInstall(install):
+    """Ensure build_ext runs first in install command."""
+    def run(self):
+        self.run_command('build_ext')
+        install.run(self)
+
+class CustomBuildExt(build_ext):
+    """ Check if openbabel bindings are installed || build them """
+    def run(self):
+        try: import openbabel
+        except ModuleNotFoundError:
+            try: import requests
+            except ModuleNotFoundError:
+                install_pkg_via_pip('requests')
+            build_ob_py_bindings()
+        return
+
 setup(name='plip',
-      version=config.__version__,
-      description='PLIP - Fully automated protein-ligand interaction profiler',
-      classifiers=[
-          'Development Status :: 5 - Production/Stable',
-          'Intended Audience :: Science/Research',
-          'Natural Language :: English',
-          'License :: OSI Approved :: GNU General Public License v2 (GPLv2)',
-          'Programming Language :: Python :: 3.6',
-          'Topic :: Scientific/Engineering :: Bio-Informatics'
-      ],
-      url='https://github.com/pharmai/plip',
-      author='PharmAI GmbH',
-      author_email='hello at pharm.ai',
-      license='GPLv2',
-      packages=find_packages(),
-      scripts=['plip/plipcmd.py'],
-      install_requires=[
-          'openbabel',
-          'numpy',
-          'lxml'
-      ],
-      entry_points={
-          "console_scripts": [
-              "plip = plip.plipcmd:main"
-          ]
-      },
-      zip_safe=False)
+    version=config.__version__,
+    description='PLIP - Fully automated protein-ligand interaction profiler',
+    classifiers=[
+        'Development Status :: 5 - Production/Stable',
+        'Intended Audience :: Science/Research',
+        'Natural Language :: English',
+        'License :: OSI Approved :: GNU General Public License v2 (GPLv2)',
+        'Programming Language :: Python :: 3.6',
+        'Topic :: Scientific/Engineering :: Bio-Informatics'
+        ],
+    url='https://github.com/pharmai/plip',
+    author='PharmAI GmbH',
+    author_email='hello at pharm.ai',
+    license='GPLv2',
+    packages=find_packages(),
+    scripts=['plip/plipcmd.py'],
+    cmdclass={'build': CustomBuild, 'build_ext': CustomBuildExt, 'install': CustomInstall},
+    install_requires=[
+        'numpy',
+        'lxml'
+        ],
+    entry_points={
+        "console_scripts": [
+            "plip = plip.plipcmd:main"
+            ]
+        },
+    zip_safe=False)



View it on GitLab: https://salsa.debian.org/med-team/plip/-/compare/9408f2de4e28772642baaea5c38da321a843cb5e...987142df85e4487ee8693e3417f93cd5d85a962c

-- 
View it on GitLab: https://salsa.debian.org/med-team/plip/-/compare/9408f2de4e28772642baaea5c38da321a843cb5e...987142df85e4487ee8693e3417f93cd5d85a962c
You're receiving this email because of your account on salsa.debian.org.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20250114/b89831d7/attachment-0001.htm>


More information about the debian-med-commit mailing list