import itertools import os import re import tempfile from collections import namedtuple from operator import itemgetter import numpy as np from openbabel import pybel from plip.basic import config, logger from plip.basic.supplemental import centroid, tilde_expansion, tmpfile, classify_by_name from plip.basic.supplemental import cluster_doubles, is_lig, normalize_vector, vector, ring_is_planar 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.structure.detection import halogen, pication, water_bridges, metal_complexation from plip.structure.detection import hydrophobic_interactions, pistacking, hbonds, saltbridge logger = logger.get_logger() class PDBParser: def __init__(self, pdbpath, as_string): self.as_string = as_string self.pdbpath = pdbpath self.num_fixed_lines = 0 self.covlinkage = namedtuple("covlinkage", "id1 chain1 pos1 conf1 id2 chain2 pos2 conf2") self.proteinmap, self.modres, self.covalent, self.altconformations, self.corrected_pdb = self.parse_pdb() def parse_pdb(self): """Extracts additional information from PDB files. I. When reading in a PDB file, OpenBabel numbers ATOMS and HETATOMS continously. In PDB files, TER records are also counted, leading to a different numbering system. This functions reads in a PDB file and provides a mapping as a dictionary. II. Additionally, it returns a list of modified residues. III. Furthermore, covalent linkages between ligands and protein residues/other ligands are identified IV. Alternative conformations """ if self.as_string: fil = self.pdbpath.rstrip('\n').split('\n') # Removing trailing newline character else: f = read(self.pdbpath) fil = f.readlines() f.close() corrected_lines = [] i, j = 0, 0 # idx and PDB numbering d = {} modres = set() covalent = [] alt = [] previous_ter = False # Standard without fixing if not config.NOFIX: if not config.PLUGIN_MODE: lastnum = 0 # Atom numbering (has to be consecutive) other_models = False for line in fil: if not other_models: # Only consider the first model in an NRM structure corrected_line, newnum = self.fix_pdbline(line, lastnum) if corrected_line is not None: if corrected_line.startswith('MODEL'): try: # Get number of MODEL (1,2,3) model_num = int(corrected_line[10:14]) if model_num > 1: # MODEL 2,3,4 etc. other_models = True except ValueError: logger.debug(f'ignoring invalid MODEL entry: {corrected_line}') corrected_lines.append(corrected_line) lastnum = newnum corrected_pdb = ''.join(corrected_lines) else: corrected_pdb = self.pdbpath corrected_lines = fil else: corrected_pdb = self.pdbpath corrected_lines = fil for line in corrected_lines: if line.startswith(("ATOM", "HETATM")): # Retrieve alternate conformations atomid, location = int(line[6:11]), line[16] location = 'A' if location == ' ' else location if location != 'A': alt.append(atomid) if not previous_ter: i += 1 j += 1 else: i += 1 j += 2 d[i] = j previous_ter = False # Numbering Changes at TER records if line.startswith("TER"): previous_ter = True # Get modified residues if line.startswith("MODRES"): modres.add(line[12:15].strip()) # Get covalent linkages between ligands if line.startswith("LINK"): covalent.append(self.get_linkage(line)) return d, modres, covalent, alt, corrected_pdb def fix_pdbline(self, pdbline, lastnum): """Fix a PDB line if information is missing.""" pdbqt_conversion = { "HD": "H", "HS": "H", "NA": "N", "NS": "N", "OA": "O", "OS": "O", "SA": "S"} fixed = False new_num = 0 forbidden_characters = "[^a-zA-Z0-9_]" pdbline = pdbline.strip('\n') # Some MD / Docking tools produce empty lines, leading to segfaults if len(pdbline.strip()) == 0: self.num_fixed_lines += 1 return None, lastnum if len(pdbline) > 100: # Should be 80 long self.num_fixed_lines += 1 return None, lastnum # TER Entries also have continuing numbering, consider them as well if pdbline.startswith('TER'): new_num = lastnum + 1 if pdbline.startswith('ATOM'): new_num = lastnum + 1 current_num = int(pdbline[6:11]) resnum = pdbline[22:27].strip() resname = pdbline[17:21].strip() # Invalid residue number try: int(resnum) except ValueError: pdbline = pdbline[:22] + ' 0 ' + pdbline[27:] fixed = True # Invalid characters in residue name if re.match(forbidden_characters, resname.strip()): pdbline = pdbline[:17] + 'UNK ' + pdbline[21:] fixed = True if lastnum + 1 != current_num: pdbline = pdbline[:6] + (5 - len(str(new_num))) * ' ' + str(new_num) + ' ' + pdbline[12:] fixed = True # No chain assigned if pdbline[21] == ' ': pdbline = pdbline[:21] + 'A' + pdbline[22:] fixed = True if pdbline.endswith('H'): self.num_fixed_lines += 1 return None, lastnum # Sometimes, converted PDB structures contain PDBQT atom types. Fix that. for pdbqttype in pdbqt_conversion: if pdbline.strip().endswith(pdbqttype): pdbline = pdbline.strip()[:-2] + ' ' + pdbqt_conversion[pdbqttype] + '\n' self.num_fixed_lines += 1 if pdbline.startswith('HETATM'): new_num = lastnum + 1 try: current_num = int(pdbline[6:11]) except ValueError: current_num = None logger.debug(f'invalid HETATM entry: {pdbline}') if lastnum + 1 != current_num: pdbline = pdbline[:6] + (5 - len(str(new_num))) * ' ' + str(new_num) + ' ' + pdbline[12:] fixed = True # No chain assigned or number assigned as chain if pdbline[21] == ' ': pdbline = pdbline[:21] + 'Z' + pdbline[22:] fixed = True # No residue number assigned if pdbline[23:26] == ' ': pdbline = pdbline[:23] + '999' + pdbline[26:] fixed = True # Non-standard Ligand Names ligname = pdbline[17:21].strip() if len(ligname) > 3: pdbline = pdbline[:17] + ligname[:3] + ' ' + pdbline[21:] fixed = True if re.match(forbidden_characters, ligname.strip()): pdbline = pdbline[:17] + 'LIG ' + pdbline[21:] fixed = True if len(ligname.strip()) == 0: pdbline = pdbline[:17] + 'LIG ' + pdbline[21:] fixed = True if pdbline.endswith('H'): self.num_fixed_lines += 1 return None, lastnum # Sometimes, converted PDB structures contain PDBQT atom types. Fix that. for pdbqttype in pdbqt_conversion: if pdbline.strip().endswith(pdbqttype): pdbline = pdbline.strip()[:-2] + ' ' + pdbqt_conversion[pdbqttype] + ' ' self.num_fixed_lines += 1 self.num_fixed_lines += 1 if fixed else 0 return pdbline + '\n', max(new_num, lastnum) def get_linkage(self, line): """Get the linkage information from a LINK entry PDB line.""" conf1, id1, chain1, pos1 = line[16].strip(), line[17:20].strip(), line[21].strip(), int(line[22:26]) conf2, id2, chain2, pos2 = line[46].strip(), line[47:50].strip(), line[51].strip(), int(line[52:56]) return self.covlinkage(id1=id1, chain1=chain1, pos1=pos1, conf1=conf1, id2=id2, chain2=chain2, pos2=pos2, conf2=conf2) class LigandFinder: def __init__(self, proteincomplex, altconf, modres, covalent, mapper): self.lignames_all = None self.lignames_kept = None self.water = None self.proteincomplex = proteincomplex self.altconformations = altconf self.modresidues = modres self.covalent = covalent self.mapper = mapper self.ligands = self.getligs() self.excluded = sorted(list(self.lignames_all.difference(set(self.lignames_kept)))) def getpeptides(self, chain): """If peptide ligand chains are defined via the command line options, try to extract the underlying ligand formed by all residues in the 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 if len(all_from_chain) == 0: return None else: non_water = [o for o in all_from_chain if not o.GetResidueProperty(9)] ligand = self.extract_ligand(non_water) return ligand def getligs(self): """Get all ligands from a PDB file and prepare them for analysis. Returns all non-empty ligands. """ if config.PEPTIDES == [] and config.INTRA is None: # Extract small molecule ligands (default) ligands = [] # Filter for ligands using lists ligand_residues, self.lignames_all, self.water = self.filter_for_ligands() all_res_dict = {(a.GetName(), a.GetChain(), a.GetNum()): a for a in ligand_residues} self.lignames_kept = list(set([a.GetName() for a in ligand_residues])) if not config.BREAKCOMPOSITE: # Update register of covalent links with those between DNA/RNA subunits self.covalent += nucleotide_linkage(all_res_dict) # Find fragment linked by covalent bonds res_kmers = self.identify_kmers(all_res_dict) else: res_kmers = [[a, ] for a in ligand_residues] logger.debug(f'{len(res_kmers)} ligand kmer(s) detected for closer inspection') for kmer in res_kmers: # iterate over all ligands and extract molecules + information if len(kmer) > config.MAX_COMPOSITE_LENGTH: logger.debug( f'ligand kmer(s) filtered out with a length of {len(kmer)} fragments ({config.MAX_COMPOSITE_LENGTH} allowed)') else: ligands.append(self.extract_ligand(kmer)) 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: peptide_ligands = [self.getpeptides(chain) for chain in config.PEPTIDES] elif config.INTRA is not None: peptide_ligands = [self.getpeptides(config.INTRA), ] ligands = [p for p in peptide_ligands if p is not None] self.covalent, self.lignames_kept, self.lignames_all = [], [], set() return [lig for lig in ligands if len(lig.mol.atoms) != 0] def extract_ligand(self, kmer): """Extract the ligand by copying atoms and bonds and assign all information necessary for later steps.""" data = namedtuple('ligand', 'mol hetid chain position water members longname type atomorder can_to_pdb') members = [(res.GetName(), res.GetChain(), int32_to_negative(res.GetNum())) for res in kmer] members = sort_members_by_importance(members) rname, rchain, rnum = members[0] logger.debug(f'finalizing extraction for ligand {rname}:{rchain}:{rnum} with {len(kmer)} elements') names = [x[0] for x in members] longname = '-'.join([x[0] for x in members]) if config.PEPTIDES: ligtype = 'PEPTIDE' elif config.INTRA is not None: ligtype = 'INTRA' else: # Classify a ligand by its HETID(s) ligtype = classify_by_name(names) logger.debug(f'ligand classified as {ligtype}') hetatoms = set() for obresidue in kmer: hetatoms_res = set([(obatom.GetIdx(), obatom) for obatom in pybel.ob.OBResidueAtomIter(obresidue) if obatom.GetAtomicNum() != 1]) if not config.ALTLOC: # Remove alternative conformations (standard -> True) hetatoms_res = set([atm for atm in hetatoms_res if not self.mapper.mapid(atm[0], mtype='protein', to='internal') in self.altconformations]) hetatoms.update(hetatoms_res) logger.debug(f'hetero atoms determined (n={len(hetatoms)})') hetatoms = dict(hetatoms) # make it a dict with idx as key and OBAtom as value lig = pybel.ob.OBMol() # new ligand mol 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') ############################################################## # 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)])) 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) lig.data.update({'Name': rname, 'Chain': rchain, 'ResNr': rnum}) # Check if a negative residue number is represented as a 32 bit integer if rnum > 10 ** 5: rnum = int32_to_negative(rnum) lig.title = ':'.join((rname, rchain, str(rnum))) self.mapper.ligandmaps[lig.title] = mapold logger.debug('renumerated molecule generated') if not config.NOPDBCANMAP: atomorder = canonicalize(lig) else: atomorder = None can_to_pdb = {} if atomorder is not None: can_to_pdb = {atomorder[key - 1]: mapold[key] for key in mapold} ligand = data(mol=lig, hetid=rname, chain=rchain, position=rnum, water=self.water, members=members, longname=longname, type=ligtype, atomorder=atomorder, can_to_pdb=can_to_pdb) return ligand def is_het_residue(self, obres): """Given an OBResidue, determines if the residue is indeed a possible ligand in the PDB file""" if not obres.GetResidueProperty(0): # If the residue is NOT amino (0) # It can be amino_nucleo, coenzme, ion, nucleo, protein, purine, pyrimidine, solvent # In these cases, it is a ligand candidate return True else: # Here, the residue is classified as amino # Amino acids can still be ligands, so we check for HETATM entries # Only residues with at least one HETATM entry are processed as ligands het_atoms = [] for atm in pybel.ob.OBResidueAtomIter(obres): het_atoms.append(obres.IsHetAtom(atm)) if True in het_atoms: return True return False def filter_for_ligands(self): """Given an OpenBabel Molecule, get all ligands, their names, and water""" candidates1 = [o for o in pybel.ob.OBResidueIter( self.proteincomplex.OBMol) if not o.GetResidueProperty(9) and self.is_het_residue(o)] if config.DNARECEPTOR: # If DNA is the receptor, don't consider DNA as a ligand candidates1 = [res for res in candidates1 if res.GetName() not in config.DNA + config.RNA] all_lignames = set([a.GetName() for a in candidates1]) water = [o for o in pybel.ob.OBResidueIter(self.proteincomplex.OBMol) if o.GetResidueProperty(9)] # Filter out non-ligands if not config.KEEPMOD: # Keep modified residues as ligands candidates2 = [a for a in candidates1 if is_lig(a.GetName()) and a.GetName() not in self.modresidues] else: candidates2 = [a for a in candidates1 if is_lig(a.GetName())] logger.debug(f'{len(candidates2)} ligand(s) after first filtering step') ############################################ # Filtering by counting and artifacts list # ############################################ artifacts = [] unique_ligs = set(a.GetName() for a in candidates2) for ulig in unique_ligs: # Discard if appearing 15 times or more and is possible artifact if ulig in config.biolip_list and [a.GetName() for a in candidates2].count(ulig) >= 15: artifacts.append(ulig) selected_ligands = [a for a in candidates2 if a.GetName() not in artifacts] return selected_ligands, all_lignames, water def identify_kmers(self, residues): """Using the covalent linkage information, find out which fragments/subunits form a ligand.""" # Remove all those not considered by ligands and pairings including alternate conformations ligdoubles = [[(link.id1, link.chain1, link.pos1), (link.id2, link.chain2, link.pos2)] for link in [c for c in self.covalent if c.id1 in self.lignames_kept and c.id2 in self.lignames_kept and c.conf1 in ['A', ''] and c.conf2 in ['A', ''] and (c.id1, c.chain1, c.pos1) in residues and (c.id2, c.chain2, c.pos2) in residues]] kmers = cluster_doubles(ligdoubles) if not kmers: # No ligand kmers, just normal independent ligands return [[residues[res]] for res in residues] else: # res_kmers contains clusters of covalently bound ligand residues (kmer ligands) res_kmers = [[residues[res] for res in kmer] for kmer in kmers] # In this case, add other ligands which are not part of a kmer in_kmer = [] for res_kmer in res_kmers: for res in res_kmer: in_kmer.append((res.GetName(), res.GetChain(), res.GetNum())) for res in residues: if res not in in_kmer: newres = [residues[res], ] res_kmers.append(newres) return res_kmers class Mapper: """Provides functions for mapping atom IDs in the correct way""" def __init__(self): self.proteinmap = None # Map internal atom IDs of protein residues to original PDB Atom IDs self.ligandmaps = {} # Map IDs of new ligand molecules to internal IDs (or PDB IDs?) self.original_structure = None def mapid(self, idx, mtype, bsid=None, to='original'): # Mapping to original IDs is standard for ligands if mtype == 'reversed': # Needed to map internal ID back to original protein ID return self.reversed_proteinmap[idx] if mtype == 'protein': return self.proteinmap[idx] elif mtype == 'ligand': if to == 'internal': return self.ligandmaps[bsid][idx] elif to == 'original': return self.proteinmap[self.ligandmaps[bsid][idx]] def id_to_atom(self, idx): """Returns the atom for a given original ligand ID. To do this, the ID is mapped to the protein first and then the atom returned. """ mapped_idx = self.mapid(idx, 'reversed') return pybel.Atom(self.original_structure.GetAtom(mapped_idx)) class Mol: def __init__(self, altconf, mapper, mtype, bsid): self.mtype = mtype self.bsid = bsid self.rings = None self.hydroph_atoms = None self.charged = None self.hbond_don_atom_pairs = None self.hbond_acc_atoms = None self.altconf = altconf self.Mapper = mapper def hydrophobic_atoms(self, all_atoms): """Select all carbon atoms which have only carbons and/or hydrogens as direct neighbors.""" atom_set = [] data = namedtuple('hydrophobic', 'atom orig_atom orig_idx') atm = [a for a in all_atoms if a.atomicnum == 6 and set([natom.GetAtomicNum() for natom in pybel.ob.OBAtomAtomIter(a.OBAtom)]).issubset( {1, 6})] for atom in atm: orig_idx = self.Mapper.mapid(atom.idx, mtype=self.mtype, bsid=self.bsid) orig_atom = self.Mapper.id_to_atom(orig_idx) if atom.idx not in self.altconf: atom_set.append(data(atom=atom, orig_atom=orig_atom, orig_idx=orig_idx)) return atom_set def find_hba(self, all_atoms): """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): 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')) a_set = sorted(a_set, key=lambda x: x.a_orig_idx) return a_set def find_hbd(self, all_atoms, hydroph_atoms): """Find all possible strong and weak hydrogen bonds donors (all hydrophobic C-H pairings)""" 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 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) d_orig_atom = self.Mapper.id_to_atom(d_orig_idx) donor_pairs.append(data(d=carbon, d_orig_atom=d_orig_atom, d_orig_idx=d_orig_idx, h=pybel.Atom(adj_atom), type='weak')) donor_pairs = sorted(donor_pairs, key=lambda x: (x.d_orig_idx, x.h.idx)) return donor_pairs def find_rings(self, mol, all_atoms): """Find rings and return only aromatic. Rings have to be sufficiently planar OR be detected by OpenBabel as aromatic.""" data = namedtuple('aromatic_ring', 'atoms orig_atoms atoms_orig_idx normal obj center type') rings = [] aromatic_amino = ['TYR', 'TRP', 'HIS', 'PHE'] ring_candidates = mol.OBMol.GetSSSR() logger.debug(f'number of aromatic ring candidates: {len(ring_candidates)}') # Check here first for ligand rings not being detected as aromatic by Babel and check for planarity for ring in ring_candidates: r_atoms = [a for a in all_atoms if ring.IsMember(a.OBAtom)] r_atoms = sorted(r_atoms, key=lambda x: x.idx) if 4 < len(r_atoms) <= 6: res = list(set([whichrestype(a) for a in r_atoms])) # re-sort ring atoms for only ligands, because HETATM numbering is not canonical in OpenBabel if res[0] == 'UNL': ligand_orig_idx = [self.Mapper.ligandmaps[self.bsid][a.idx] for a in r_atoms] sort_order = np.argsort(np.array(ligand_orig_idx)) r_atoms = [r_atoms[i] for i in sort_order] if ring.IsAromatic() or res[0] in aromatic_amino or ring_is_planar(ring, r_atoms): # Causes segfault with OpenBabel 2.3.2, so deactivated # typ = ring.GetType() if not ring.GetType() == '' else 'unknown' # Alternative typing ring_type = '%s-membered' % len(r_atoms) ring_atms = [r_atoms[a].coords for a in [0, 2, 4]] # Probe atoms for normals, assuming planarity ringv1 = vector(ring_atms[0], ring_atms[1]) ringv2 = vector(ring_atms[2], ring_atms[0]) atoms_orig_idx = [self.Mapper.mapid(r_atom.idx, mtype=self.mtype, bsid=self.bsid) for r_atom in r_atoms] orig_atoms = [self.Mapper.id_to_atom(idx) for idx in atoms_orig_idx] rings.append(data(atoms=r_atoms, orig_atoms=orig_atoms, atoms_orig_idx=atoms_orig_idx, normal=normalize_vector(np.cross(ringv1, ringv2)), obj=ring, center=centroid([ra.coords for ra in r_atoms]), type=ring_type)) return rings def get_hydrophobic_atoms(self): return self.hydroph_atoms def get_hba(self): return self.hbond_acc_atoms def get_hbd(self): return [don_pair for don_pair in self.hbond_don_atom_pairs if don_pair.type == 'regular'] def get_weak_hbd(self): return [don_pair for don_pair in self.hbond_don_atom_pairs if don_pair.type == 'weak'] def get_pos_charged(self): return [charge for charge in self.charged if charge.type == 'positive'] def get_neg_charged(self): return [charge for charge in self.charged if charge.type == 'negative'] class PLInteraction: """Class to store a ligand, a protein and their interactions.""" def __init__(self, lig_obj, bs_obj, protcomplex): """Detect all interactions when initializing""" self.ligand = lig_obj self.lig_members = lig_obj.members self.pdbid = protcomplex.pymol_name self.bindingsite = bs_obj self.Mapper = protcomplex.Mapper self.output_path = protcomplex.output_path self.altconf = protcomplex.altconf # #@todo Refactor code to combine different directionality self.saltbridge_lneg = saltbridge(self.bindingsite.get_pos_charged(), self.ligand.get_neg_charged(), True) self.saltbridge_pneg = saltbridge(self.ligand.get_pos_charged(), self.bindingsite.get_neg_charged(), False) self.all_hbonds_ldon = hbonds(self.bindingsite.get_hba(), self.ligand.get_hbd(), False, 'strong') self.all_hbonds_pdon = hbonds(self.ligand.get_hba(), self.bindingsite.get_hbd(), True, 'strong') self.hbonds_ldon = self.refine_hbonds_ldon(self.all_hbonds_ldon, self.saltbridge_lneg, self.saltbridge_pneg) self.hbonds_pdon = self.refine_hbonds_pdon(self.all_hbonds_pdon, self.saltbridge_lneg, self.saltbridge_pneg) 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.pication_laro = self.refine_pi_cation_laro(self.all_pi_cation_laro, self.pistacking) self.all_hydrophobic_contacts = hydrophobic_interactions(self.bindingsite.get_hydrophobic_atoms(), self.ligand.get_hydrophobic_atoms()) self.hydrophobic_contacts = self.refine_hydrophobic(self.all_hydrophobic_contacts, self.pistacking) self.halogen_bonds = halogen(self.bindingsite.halogenbond_acc, self.ligand.halogenbond_don) self.water_bridges = water_bridges(self.bindingsite.get_hba(), self.ligand.get_hba(), self.bindingsite.get_hbd(), self.ligand.get_hbd(), self.ligand.water) self.water_bridges = self.refine_water_bridges(self.water_bridges, self.hbonds_ldon, self.hbonds_pdon) self.metal_complexes = metal_complexation(self.ligand.metals, self.ligand.metal_binding, self.bindingsite.metal_binding) self.all_itypes = self.saltbridge_lneg + self.saltbridge_pneg + self.hbonds_pdon self.all_itypes = self.all_itypes + self.hbonds_ldon + self.pistacking + self.pication_laro + self.pication_paro self.all_itypes = self.all_itypes + self.hydrophobic_contacts + self.halogen_bonds + self.water_bridges self.all_itypes = self.all_itypes + self.metal_complexes self.no_interactions = all(len(i) == 0 for i in self.all_itypes) self.unpaired_hba, self.unpaired_hbd, self.unpaired_hal = self.find_unpaired_ligand() self.unpaired_hba_orig_idx = [self.Mapper.mapid(atom.idx, mtype='ligand', bsid=self.ligand.bsid) for atom in self.unpaired_hba] self.unpaired_hbd_orig_idx = [self.Mapper.mapid(atom.idx, mtype='ligand', bsid=self.ligand.bsid) for atom in self.unpaired_hbd] self.unpaired_hal_orig_idx = [self.Mapper.mapid(atom.idx, mtype='ligand', bsid=self.ligand.bsid) for atom in self.unpaired_hal] self.num_unpaired_hba, self.num_unpaired_hbd = len(self.unpaired_hba), len(self.unpaired_hbd) self.num_unpaired_hal = len(self.unpaired_hal) # Exclude empty chains (coming from ligand as a target, from metal complexes) self.interacting_chains = sorted(list(set([i.reschain for i in self.all_itypes if i.reschain not in [' ', None]]))) # Get all interacting residues, excluding ligand and water molecules self.interacting_res = list(set([''.join([str(i.resnr), str(i.reschain)]) for i in self.all_itypes if i.restype not in ['LIG', 'HOH']])) if len(self.interacting_res) != 0: logger.info( f'ligand interacts with {len(self.interacting_res)} binding site residue(s) in chain(s) {self.interacting_chains}') interactions_list = [] num_saltbridges = len(self.saltbridge_lneg + self.saltbridge_pneg) num_hbonds = len(self.hbonds_ldon + self.hbonds_pdon) num_pication = len(self.pication_laro + self.pication_paro) num_pistack = len(self.pistacking) num_halogen = len(self.halogen_bonds) num_waterbridges = len(self.water_bridges) if num_saltbridges != 0: interactions_list.append('%i salt bridge(s)' % num_saltbridges) if num_hbonds != 0: interactions_list.append('%i hydrogen bond(s)' % num_hbonds) if num_pication != 0: interactions_list.append('%i pi-cation interaction(s)' % num_pication) if num_pistack != 0: interactions_list.append('%i pi-stacking(s)' % num_pistack) if num_halogen != 0: interactions_list.append('%i halogen bond(s)' % num_halogen) if num_waterbridges != 0: interactions_list.append('%i water bridge(s)' % num_waterbridges) if not len(interactions_list) == 0: logger.info(f'complex uses {interactions_list}') else: logger.info('no interactions for this ligand') def find_unpaired_ligand(self): """Identify unpaired functional in groups in ligands, involving H-Bond donors, acceptors, halogen bond donors. """ unpaired_hba, unpaired_hbd, unpaired_hal = [], [], [] # Unpaired hydrogen bond acceptors/donors in ligand (not used for hydrogen bonds/water, salt bridges/mcomplex) involved_atoms = [hbond.a.idx for hbond in self.hbonds_pdon] + [hbond.d.idx for hbond in self.hbonds_ldon] [[involved_atoms.append(atom.idx) for atom in sb.negative.atoms] for sb in self.saltbridge_lneg] [[involved_atoms.append(atom.idx) for atom in sb.positive.atoms] for sb in self.saltbridge_pneg] [involved_atoms.append(wb.a.idx) for wb in self.water_bridges if wb.protisdon] [involved_atoms.append(wb.d.idx) for wb in self.water_bridges if not wb.protisdon] [involved_atoms.append(mcomplex.target.atom.idx) for mcomplex in self.metal_complexes if mcomplex.location == 'ligand'] for atom in [hba.a for hba in self.ligand.get_hba()]: if atom.idx not in involved_atoms: unpaired_hba.append(atom) for atom in [hbd.d for hbd in self.ligand.get_hbd()]: if atom.idx not in involved_atoms: unpaired_hbd.append(atom) # unpaired halogen bond donors in ligand (not used for the previous + halogen bonds) [involved_atoms.append(atom.don.x.idx) for atom in self.halogen_bonds] for atom in [haldon.x for haldon in self.ligand.halogenbond_don]: if atom.idx not in involved_atoms: unpaired_hal.append(atom) return unpaired_hba, unpaired_hbd, unpaired_hal def refine_hydrophobic(self, all_h, pistacks): """Apply several rules to reduce the number of hydrophobic interactions.""" sel = {} # 1. Rings interacting via stacking can't have additional hydrophobic contacts between each other. for pistack, h in itertools.product(pistacks, all_h): h1, h2 = h.bsatom.idx, h.ligatom.idx brs, lrs = [p1.idx for p1 in pistack.proteinring.atoms], [p2.idx for p2 in pistack.ligandring.atoms] if h1 in brs and h2 in lrs: sel[(h1, h2)] = "EXCLUDE" hydroph = [h for h in all_h if not (h.bsatom.idx, h.ligatom.idx) in sel] sel2 = {} # 2. If a ligand atom interacts with several binding site atoms in the same residue, # keep only the one with the closest distance for h in hydroph: if not (h.ligatom.idx, h.resnr) in sel2: sel2[(h.ligatom.idx, h.resnr)] = h else: if sel2[(h.ligatom.idx, h.resnr)].distance > h.distance: 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) 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}') return hydroph_final def refine_hbonds_ldon(self, all_hbonds, salt_lneg, salt_pneg): """Refine selection of hydrogen bonds. Do not allow groups which already form salt bridges to form H-Bonds.""" i_set = {} for hbond in all_hbonds: i_set[hbond] = False for salt in salt_pneg: protidx, ligidx = [at.idx for at in salt.negative.atoms], [at.idx for at in salt.positive.atoms] if hbond.d.idx in ligidx and hbond.a.idx in protidx: i_set[hbond] = True for salt in salt_lneg: protidx, ligidx = [at.idx for at in salt.positive.atoms], [at.idx for at in salt.negative.atoms] if hbond.d.idx in ligidx and hbond.a.idx in protidx: i_set[hbond] = True # Allow only one hydrogen bond per donor, select interaction with larger donor angle second_set = {} hbls = [k for k in i_set.keys() if not i_set[k]] for hbl in hbls: if hbl.d.idx not in second_set: second_set[hbl.d.idx] = (hbl.angle, hbl) else: if second_set[hbl.d.idx][0] < hbl.angle: second_set[hbl.d.idx] = (hbl.angle, hbl) return [hb[1] for hb in second_set.values()] def refine_hbonds_pdon(self, all_hbonds, salt_lneg, salt_pneg): """Refine selection of hydrogen bonds. Do not allow groups which already form salt bridges to form H-Bonds with atoms of the same group. """ i_set = {} for hbond in all_hbonds: i_set[hbond] = False for salt in salt_lneg: protidx, ligidx = [at.idx for at in salt.positive.atoms], [at.idx for at in salt.negative.atoms] if hbond.a.idx in ligidx and hbond.d.idx in protidx: i_set[hbond] = True for salt in salt_pneg: protidx, ligidx = [at.idx for at in salt.negative.atoms], [at.idx for at in salt.positive.atoms] if hbond.a.idx in ligidx and hbond.d.idx in protidx: i_set[hbond] = True # Allow only one hydrogen bond per donor, select interaction with larger donor angle second_set = {} hbps = [k for k in i_set.keys() if not i_set[k]] for hbp in hbps: if hbp.d.idx not in second_set: second_set[hbp.d.idx] = (hbp.angle, hbp) else: if second_set[hbp.d.idx][0] < hbp.angle: second_set[hbp.d.idx] = (hbp.angle, hbp) return [hb[1] for hb in second_set.values()] def refine_pi_cation_laro(self, 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. """ 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 not exclude: i_set.append(picat) return i_set def refine_water_bridges(self, 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] 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 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 for wb_tuple in wb_dict: water, acceptor = wb_tuple if water not in wb_dict2: wb_dict2[water] = [(abs(omega - wb_dict[wb_tuple].w_angle), wb_dict[wb_tuple]), ] elif len(wb_dict2[water]) == 1: 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): wb_dict2[water] = [wb_dict2[water][0], (wb_dict[wb_tuple].w_angle, wb_dict[wb_tuple])] filtered_wb = [] for fwbridges in wb_dict2.values(): [filtered_wb.append(fwb[1]) for fwb in fwbridges] return filtered_wb class BindingSite(Mol): def __init__(self, atoms, protcomplex, cclass, altconf, min_dist, mapper): """Find all relevant parts which could take part in interactions""" Mol.__init__(self, altconf, mapper, mtype='protein', bsid=None) self.complex = cclass self.full_mol = protcomplex self.all_atoms = atoms self.min_dist = min_dist # Minimum distance of bs res to ligand self.bs_res = list(set([''.join([str(whichresnumber(a)), whichchain(a)]) for a in self.all_atoms])) # e.g. 47A self.rings = self.find_rings(self.full_mol, self.all_atoms) self.hydroph_atoms = self.hydrophobic_atoms(self.all_atoms) self.hbond_acc_atoms = self.find_hba(self.all_atoms) self.hbond_don_atom_pairs = self.find_hbd(self.all_atoms, self.hydroph_atoms) self.charged = self.find_charged(self.full_mol) self.halogenbond_acc = self.find_hal(self.all_atoms) self.metal_binding = self.find_metal_binding(self.full_mol) def find_hal(self, atoms): """Look for halogen bond acceptors (Y-{O|P|N|S}, with Y=C,P,S)""" data = namedtuple('hal_acceptor', 'o o_orig_idx y y_orig_idx') a_set = [] # All oxygens, nitrogen, sulfurs with neighboring carbon, phosphor, nitrogen or sulfur for a in [at for at in atoms if at.atomicnum in [8, 7, 16]]: n_atoms = [na for na in pybel.ob.OBAtomAtomIter(a.OBAtom) if na.GetAtomicNum() in [6, 7, 15, 16]] if len(n_atoms) == 1: # Proximal atom o_orig_idx = self.Mapper.mapid(a.idx, mtype=self.mtype, bsid=self.bsid) y_orig_idx = self.Mapper.mapid(n_atoms[0].GetIdx(), mtype=self.mtype, bsid=self.bsid) a_set.append(data(o=a, o_orig_idx=o_orig_idx, y=pybel.Atom(n_atoms[0]), y_orig_idx=y_orig_idx)) return a_set def find_charged(self, mol): """Looks for positive charges in arginine, histidine or lysine, for negative in aspartic and glutamic acid.""" 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]: 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='protein') in self.altconf: 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='positive', center=centroid([ac.coords for ac in a_contributing]), restype=res.GetName(), resnr=res.GetNum(), reschain=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='protein') in self.altconf: 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', center=centroid([ac.coords for ac in a_contributing]), restype=res.GetName(), resnr=res.GetNum(), reschain=res.GetChain())) return a_set def find_metal_binding(self, mol): """Looks for atoms that could possibly be involved in chelating a metal ion. This can be any main chain oxygen atom or oxygen, nitrogen and sulfur from specific amino acids""" data = namedtuple('metal_binding', 'atom atom_orig_idx type restype resnr reschain location') a_set = [] for res in pybel.ob.OBResidueIter(mol.OBMol): restype, reschain, resnr = res.GetName().upper(), res.GetChain(), res.GetNum() if restype in ['ASP', 'GLU', 'SER', 'THR', 'TYR']: # Look for oxygens here 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='protein') in self.altconf: atom_orig_idx = self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid) a_set.append(data(atom=pybel.Atom(a), atom_orig_idx=atom_orig_idx, type='O', restype=restype, resnr=resnr, reschain=reschain, location='protein.sidechain')) if restype == 'HIS': # Look for nitrogen here 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='protein') in self.altconf: atom_orig_idx = self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid) a_set.append(data(atom=pybel.Atom(a), atom_orig_idx=atom_orig_idx, type='N', restype=restype, resnr=resnr, reschain=reschain, location='protein.sidechain')) if restype == 'CYS': # Look for sulfur here for a in pybel.ob.OBResidueAtomIter(res): if a.GetType().startswith('S') and res.GetAtomProperty(a, 8) \ and not self.Mapper.mapid(a.GetIdx(), mtype='protein') in self.altconf: atom_orig_idx = self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid) a_set.append(data(atom=pybel.Atom(a), atom_orig_idx=atom_orig_idx, type='S', restype=restype, resnr=resnr, reschain=reschain, location='protein.sidechain')) for a in pybel.ob.OBResidueAtomIter(res): # All main chain oxygens if a.GetType().startswith('O') and res.GetAtomProperty(a, 2) \ and not self.Mapper.mapid(a.GetIdx(), mtype='protein') in self.altconf and restype != 'HOH': atom_orig_idx = self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid) a_set.append(data(atom=pybel.Atom(a), atom_orig_idx=atom_orig_idx, type='O', restype=res.GetName(), resnr=res.GetNum(), reschain=res.GetChain(), location='protein.mainchain')) return a_set class Ligand(Mol): def __init__(self, cclass, ligand): altconf = cclass.altconf self.hetid, self.chain, self.position = ligand.hetid, ligand.chain, ligand.position self.bsid = ':'.join([self.hetid, self.chain, str(self.position)]) Mol.__init__(self, altconf, cclass.Mapper, mtype='ligand', bsid=self.bsid) self.members = ligand.members self.longname = ligand.longname self.type = ligand.type self.complex = cclass self.molecule = ligand.mol # Pybel Molecule self.smiles = self.molecule.write(format='can') # SMILES String self.inchikey = self.molecule.write(format='inchikey') self.can_to_pdb = ligand.can_to_pdb if not len(self.smiles) == 0: self.smiles = self.smiles.split()[0] else: logger.warning(f'could not write SMILES for ligand {ligand}') self.smiles = '' self.heavy_atoms = self.molecule.OBMol.NumHvyAtoms() # Heavy atoms count self.all_atoms = self.molecule.atoms self.atmdict = {l.idx: l for l in self.all_atoms} self.rings = self.find_rings(self.molecule, self.all_atoms) self.hydroph_atoms = self.hydrophobic_atoms(self.all_atoms) self.hbond_acc_atoms = self.find_hba(self.all_atoms) self.num_rings = len(self.rings) if self.num_rings != 0: logger.info(f'contains {self.num_rings} aromatic ring(s)') descvalues = self.molecule.calcdesc() self.molweight, self.logp = float(descvalues['MW']), float(descvalues['logP']) self.num_rot_bonds = int(self.molecule.OBMol.NumRotors()) self.atomorder = ligand.atomorder ########################################################## # Special Case for hydrogen bond acceptor identification # ########################################################## self.inverse_mapping = {v: k for k, v in self.Mapper.ligandmaps[self.bsid].items()} self.pdb_to_idx_mapping = {v: k for k, v in self.Mapper.proteinmap.items()} self.hbond_don_atom_pairs = self.find_hbd(self.all_atoms, self.hydroph_atoms) ###### donor_pairs = [] data = namedtuple('hbonddonor', 'd d_orig_atom d_orig_idx h type') for donor in self.all_atoms: pdbidx = self.Mapper.mapid(donor.idx, mtype='ligand', bsid=self.bsid, to='original') d = cclass.atoms[self.pdb_to_idx_mapping[pdbidx]] if d.OBAtom.IsHbondDonor(): for adj_atom in [a for a in pybel.ob.OBAtomAtomIter(d.OBAtom) if a.IsHbondDonorH()]: d_orig_atom = self.Mapper.id_to_atom(pdbidx) donor_pairs.append(data(d=donor, d_orig_atom=d_orig_atom, d_orig_idx=pdbidx, h=pybel.Atom(adj_atom), type='regular')) self.hbond_don_atom_pairs = donor_pairs ####### self.charged = self.find_charged(self.all_atoms) self.centroid = centroid([a.coords for a in self.all_atoms]) self.max_dist_to_center = max((euclidean3d(self.centroid, a.coords) for a in self.all_atoms)) self.water = [] data = namedtuple('water', 'oxy oxy_orig_idx') for hoh in ligand.water: oxy = None for at in pybel.ob.OBResidueAtomIter(hoh): if at.GetAtomicNum() == 8 and at.GetIdx() not in self.altconf: oxy = pybel.Atom(at) # There are some cases where there is no oxygen in a water residue, ignore those if not set([at.GetAtomicNum() for at in pybel.ob.OBResidueAtomIter(hoh)]) == {1} and oxy is not None: if euclidean3d(self.centroid, oxy.coords) < self.max_dist_to_center + config.BS_DIST: oxy_orig_idx = self.Mapper.mapid(oxy.idx, mtype='protein') self.water.append(data(oxy=oxy, oxy_orig_idx=oxy_orig_idx)) self.halogenbond_don = self.find_hal(self.all_atoms) self.metal_binding = self.find_metal_binding(self.all_atoms, self.water) self.metals = [] data = namedtuple('metal', 'm orig_m m_orig_idx') for a in [a for a in self.all_atoms if a.type.upper() in config.METAL_IONS]: m_orig_idx = self.Mapper.mapid(a.idx, mtype=self.mtype, bsid=self.bsid) orig_m = self.Mapper.id_to_atom(m_orig_idx) self.metals.append(data(m=a, m_orig_idx=m_orig_idx, orig_m=orig_m)) self.num_hba, self.num_hbd = len(self.hbond_acc_atoms), len(self.hbond_don_atom_pairs) self.num_hal = len(self.halogenbond_don) def get_canonical_num(self, atomnum): """Converts internal atom ID into canonical atom ID. Agrees with Canonical SMILES in XML.""" return self.atomorder[atomnum - 1] def is_functional_group(self, atom, group): """Given a pybel atom, look up if it belongs to a function group""" n_atoms = [a_neighbor.GetAtomicNum() for a_neighbor in pybel.ob.OBAtomAtomIter(atom.OBAtom)] if group in ['quartamine', 'tertamine'] and atom.atomicnum == 7: # Nitrogen # It's a nitrogen, so could be a protonated amine or quaternary ammonium if '1' not in n_atoms and len(n_atoms) == 4: return True if group == 'quartamine' else False # It's a quat. ammonium (N with 4 residues != H) elif atom.OBAtom.GetHyb() == 3 and len(n_atoms) >= 3: return True if group == 'tertamine' else False # It's sp3-hybridized, so could pick up an hydrogen else: return False if group in ['sulfonium', 'sulfonicacid', 'sulfate'] and atom.atomicnum == 16: # Sulfur if '1' not in n_atoms and len(n_atoms) == 3: # It's a sulfonium (S with 3 residues != H) return True if group == 'sulfonium' else False elif n_atoms.count(8) == 3: # It's a sulfonate or sulfonic acid return True if group == 'sulfonicacid' else False elif n_atoms.count(8) == 4: # It's a sulfate return True if group == 'sulfate' else False if group == 'phosphate' and atom.atomicnum == 15: # Phosphor if set(n_atoms) == {8}: # It's a phosphate return True if group in ['carboxylate', 'guanidine'] and atom.atomicnum == 6: # It's a carbon atom if n_atoms.count(8) == 2 and n_atoms.count(6) == 1: # It's a carboxylate group return True if group == 'carboxylate' else False elif n_atoms.count(7) == 3 and len(n_atoms) == 3: # It's a guanidine group nitro_partners = [] for nitro in pybel.ob.OBAtomAtomIter(atom.OBAtom): nitro_partners.append(len([b_neighbor for b_neighbor in pybel.ob.OBAtomAtomIter(nitro)])) if min(nitro_partners) == 1: # One nitrogen is only connected to the carbon, can pick up a H return True if group == 'guanidine' else False if group == 'halocarbon' and atom.atomicnum in [9, 17, 35, 53]: # Halogen atoms n_atoms = [na for na in pybel.ob.OBAtomAtomIter(atom.OBAtom) if na.GetAtomicNum() == 6] if len(n_atoms) == 1: # Halocarbon return True else: return False def find_hal(self, atoms): """Look for halogen bond donors (X-C, with X=F, Cl, Br, I)""" data = namedtuple('hal_donor', 'x orig_x x_orig_idx c c_orig_idx') a_set = [] for a in atoms: if self.is_functional_group(a, 'halocarbon'): n_atoms = [na for na in pybel.ob.OBAtomAtomIter(a.OBAtom) if na.GetAtomicNum() == 6] x_orig_idx = self.Mapper.mapid(a.idx, mtype=self.mtype, bsid=self.bsid) orig_x = self.Mapper.id_to_atom(x_orig_idx) c_orig_idx = [self.Mapper.mapid(na.GetIdx(), mtype=self.mtype, bsid=self.bsid) for na in n_atoms] a_set.append(data(x=a, orig_x=orig_x, x_orig_idx=x_orig_idx, c=pybel.Atom(n_atoms[0]), c_orig_idx=c_orig_idx)) if len(a_set) != 0: logger.info(f'ligand contains {len(a_set)} halogen atom(s)') return a_set def find_charged(self, all_atoms): """Identify all positively charged groups in a ligand. This search is not exhaustive, as the cases can be quite diverse. The typical cases seem to be protonated amines, quaternary ammoinium and sulfonium as mentioned in 'Cation-pi interactions in ligand recognition and catalysis' (Zacharias et al., 2002)). Identify negatively charged groups in the ligand. """ 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')) return a_set def find_metal_binding(self, lig_atoms, water_oxygens): """Looks for atoms that could possibly be involved in binding a metal ion. This can be any water oxygen, as well as oxygen from carboxylate, phophoryl, phenolate, alcohol; nitrogen from imidazole; sulfur from thiolate. """ a_set = [] data = namedtuple('metal_binding', 'atom orig_atom atom_orig_idx type fgroup restype resnr reschain location') for oxygen in water_oxygens: a_set.append(data(atom=oxygen.oxy, atom_orig_idx=oxygen.oxy_orig_idx, type='O', fgroup='water', restype=whichrestype(oxygen.oxy), resnr=whichresnumber(oxygen.oxy), reschain=whichchain(oxygen.oxy), location='water', orig_atom=self.Mapper.id_to_atom(oxygen.oxy_orig_idx))) # #@todo Refactor code for a in lig_atoms: a_orig_idx = self.Mapper.mapid(a.idx, mtype='ligand', bsid=self.bsid) n_atoms = pybel.ob.OBAtomAtomIter(a.OBAtom) # Neighboring atoms # All atomic numbers of neighboring atoms n_atoms_atomicnum = [n.GetAtomicNum() for n in pybel.ob.OBAtomAtomIter(a.OBAtom)] if a.atomicnum == 8: # Oxygen if n_atoms_atomicnum.count('1') == 1 and len(n_atoms_atomicnum) == 2: # Oxygen in alcohol (R-[O]-H) a_set.append(data(atom=a, atom_orig_idx=a_orig_idx, type='O', fgroup='alcohol', restype=self.hetid, resnr=self.position, reschain=self.chain, location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx))) if True in [n.IsAromatic() for n in n_atoms] and not a.OBAtom.IsAromatic(): # Phenolate oxygen a_set.append(data(atom=a, atom_orig_idx=a_orig_idx, type='O', fgroup='phenolate', restype=self.hetid, resnr=self.position, reschain=self.chain, location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx))) if a.atomicnum == 6: # It's a carbon atom if n_atoms_atomicnum.count(8) == 2 and n_atoms_atomicnum.count(6) == 1: # It's a carboxylate group for neighbor in [n for n in n_atoms if n.GetAtomicNum() == 8]: neighbor_orig_idx = self.Mapper.mapid(neighbor.GetIdx(), mtype='ligand', bsid=self.bsid) a_set.append(data(atom=pybel.Atom(neighbor), atom_orig_idx=neighbor_orig_idx, type='O', fgroup='carboxylate', restype=self.hetid, resnr=self.position, reschain=self.chain, location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx))) if a.atomicnum == 15: # It's a phosphor atom if n_atoms_atomicnum.count(8) >= 3: # It's a phosphoryl for neighbor in [n for n in n_atoms if n.GetAtomicNum() == 8]: neighbor_orig_idx = self.Mapper.mapid(neighbor.GetIdx(), mtype='ligand', bsid=self.bsid) a_set.append(data(atom=pybel.Atom(neighbor), atom_orig_idx=neighbor_orig_idx, type='O', fgroup='phosphoryl', restype=self.hetid, resnr=self.position, reschain=self.chain, location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx))) if n_atoms_atomicnum.count(8) == 2: # It's another phosphor-containing group #@todo (correct name?) for neighbor in [n for n in n_atoms if n.GetAtomicNum() == 8]: neighbor_orig_idx = self.Mapper.mapid(neighbor.GetIdx(), mtype='ligand', bsid=self.bsid) a_set.append(data(atom=pybel.Atom(neighbor), atom_orig_idx=neighbor_orig_idx, type='O', fgroup='phosphor.other', restype=self.hetid, resnr=self.position, reschain=self.chain, location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx))) if a.atomicnum == 7: # It's a nitrogen atom if n_atoms_atomicnum.count(6) == 2: # It's imidazole/pyrrole or similar a_set.append(data(atom=a, atom_orig_idx=a_orig_idx, type='N', fgroup='imidazole/pyrrole', restype=self.hetid, resnr=self.position, reschain=self.chain, location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx))) if a.atomicnum == 16: # It's a sulfur atom if True in [n.IsAromatic() for n in n_atoms] and not a.OBAtom.IsAromatic(): # Thiolate a_set.append(data(atom=a, atom_orig_idx=a_orig_idx, type='S', fgroup='thiolate', restype=self.hetid, resnr=self.position, reschain=self.chain, location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx))) if set(n_atoms_atomicnum) == {26}: # Sulfur in Iron sulfur cluster a_set.append(data(atom=a, atom_orig_idx=a_orig_idx, type='S', fgroup='iron-sulfur.cluster', restype=self.hetid, resnr=self.position, reschain=self.chain, location='ligand', orig_atom=self.Mapper.id_to_atom(a_orig_idx))) return a_set class PDBComplex: """Contains a collection of objects associated with a PDB complex, i.e. one or several ligands and their binding sites as well as information about the pliprofiler between them. Provides functions to load and prepare input files such as PDB files. """ def __init__(self): self.interaction_sets = {} # Dictionary with site identifiers as keys and object as value self.protcomplex = None self.filetype = None self.atoms = {} # Dictionary of Pybel atoms, accessible by their idx self.sourcefiles = {} self.information = {} self.corrected_pdb = '' self._output_path = tempfile.gettempdir() self.pymol_name = None self.modres = set() self.resis = [] self.altconf = [] # Atom idx of atoms with alternate conformations self.covalent = [] # Covalent linkages between ligands and protein residues/other ligands self.excluded = [] # Excluded ligands self.Mapper = Mapper() self.ligands = [] def __str__(self): formatted_lig_names = [":".join([x.hetid, x.chain, str(x.position)]) for x in self.ligands] return "Protein structure %s with ligands:\n" % (self.pymol_name) + "\n".join( [lig for lig in formatted_lig_names]) def load_pdb(self, pdbpath, as_string=False): """Loads a pdb file with protein AND ligand(s), separates and prepares them. If specified 'as_string', the input is a PDB string instead of a path.""" if as_string: self.sourcefiles['pdbcomplex.original'] = None self.sourcefiles['pdbcomplex'] = None self.sourcefiles['pdbstring'] = pdbpath else: self.sourcefiles['pdbcomplex.original'] = pdbpath self.sourcefiles['pdbcomplex'] = pdbpath self.information['pdbfixes'] = False pdbparser = PDBParser(pdbpath, as_string=as_string) # Parse PDB file to find errors and get additional data # #@todo Refactor and rename here self.Mapper.proteinmap = pdbparser.proteinmap self.Mapper.reversed_proteinmap = {v: k for k, v in self.Mapper.proteinmap.items()} self.modres = pdbparser.modres self.covalent = pdbparser.covalent self.altconf = pdbparser.altconformations self.corrected_pdb = pdbparser.corrected_pdb if not config.PLUGIN_MODE: if pdbparser.num_fixed_lines > 0: logger.info(f'{pdbparser.num_fixed_lines} lines automatically fixed in PDB input file') # Save modified PDB file if not as_string: basename = os.path.basename(pdbpath).split('.')[0] else: basename = "from_stdin" pdbpath_fixed = tmpfile(prefix='plipfixed.' + basename + '_', direc=self.output_path) create_folder_if_not_exists(self.output_path) self.sourcefiles['pdbcomplex'] = pdbpath_fixed self.corrected_pdb = re.sub(r'[^\x00-\x7F]+', ' ', self.corrected_pdb) # Strip non-unicode chars if not config.NOFIXFILE: # Only write to file if this option is not activated with open(pdbpath_fixed, 'w') as f: f.write(self.corrected_pdb) self.information['pdbfixes'] = True if not as_string: self.sourcefiles['filename'] = os.path.basename(self.sourcefiles['pdbcomplex']) self.protcomplex, self.filetype = read_pdb(self.corrected_pdb, as_string=True) # Update the model in the Mapper class instance self.Mapper.original_structure = self.protcomplex.OBMol logger.info('PDB structure successfully read') # Determine (temporary) PyMOL Name from Filename self.pymol_name = pdbpath.split('/')[-1].split('.')[0] + '-Protein' # Replace characters causing problems in PyMOL self.pymol_name = self.pymol_name.replace(' ', '').replace('(', '').replace(')', '').replace('-', '_') # But if possible, name it after PDBID in Header if 'HEADER' in self.protcomplex.data: # If the PDB file has a proper header potential_name = self.protcomplex.data['HEADER'][56:60].lower() if extract_pdbid(potential_name) != 'UnknownProtein': self.pymol_name = potential_name logger.debug(f'PyMOL name set as: {self.pymol_name}') # Extract and prepare ligands ligandfinder = LigandFinder(self.protcomplex, self.altconf, self.modres, self.covalent, self.Mapper) self.ligands = ligandfinder.ligands self.excluded = ligandfinder.excluded # decide whether to add polar hydrogens if not config.NOHYDRO: if not as_string: basename = os.path.basename(pdbpath).split('.')[0] else: basename = "from_stdin" self.protcomplex.OBMol.AddPolarHydrogens() output_path = os.path.join(self._output_path, f'{basename}_protonated.pdb') self.protcomplex.write('pdb', output_path, overwrite=True) logger.info(f'protonated structure written to {output_path}') else: logger.warning('no polar hydrogens will be assigned (make sure your structure contains hydrogens)') for atm in self.protcomplex: self.atoms[atm.idx] = atm if len(self.excluded) != 0: logger.info(f'excluded molecules as ligands: {self.excluded}') if config.DNARECEPTOR: self.resis = [obres for obres in pybel.ob.OBResidueIter( self.protcomplex.OBMol) if obres.GetName() in config.DNA + config.RNA] else: self.resis = [obres for obres in pybel.ob.OBResidueIter( self.protcomplex.OBMol) if obres.GetResidueProperty(0)] num_ligs = len(self.ligands) if num_ligs == 1: logger.info('analyzing one ligand') elif num_ligs > 1: logger.info(f'analyzing {num_ligs} ligands') else: logger.info(f'structure contains no ligands') def analyze(self): """Triggers analysis of all complexes in structure""" for ligand in self.ligands: self.characterize_complex(ligand) def characterize_complex(self, ligand): """Handles all basic functions for characterizing the interactions for one ligand""" single_sites = [] for member in ligand.members: single_sites.append(':'.join([str(x) for x in member])) site = ' + '.join(single_sites) site = site if not len(site) > 20 else site[:20] + '...' longname = ligand.longname if not len(ligand.longname) > 20 else ligand.longname[:20] + '...' ligtype = 'unspecified type' if ligand.type == 'UNSPECIFIED' else ligand.type ligtext = f'{longname} [{ligtype}] -- {site}' logger.info(f'processing ligand {ligtext}') if ligtype == 'PEPTIDE': logger.info(f'chain {ligand.chain} will be processed in [PEPTIDE / INTER-CHAIN] mode') if ligtype == 'INTRA': logger.info(f'chain {ligand.chain} will be processed in [INTRA-CHAIN] mode') any_in_biolip = len(set([x[0] for x in ligand.members]).intersection(config.biolip_list)) != 0 if ligtype not in ['POLYMER', 'DNA', 'ION', 'DNA+ION', 'RNA+ION', 'SMALLMOLECULE+ION'] and any_in_biolip: logger.info('may be biologically irrelevant') lig_obj = Ligand(self, ligand) cutoff = lig_obj.max_dist_to_center + config.BS_DIST bs_res = self.extract_bs(cutoff, lig_obj.centroid, self.resis) # Get a list of all atoms belonging to the binding site, search by idx bs_atoms = [self.atoms[idx] for idx in [i for i in self.atoms.keys() if self.atoms[i].OBAtom.GetResidue().GetIdx() in bs_res] if idx in self.Mapper.proteinmap and self.Mapper.mapid(idx, mtype='protein') not in self.altconf] if ligand.type == 'PEPTIDE': # If peptide, don't consider the peptide chain as part of the protein binding site bs_atoms = [a for a in bs_atoms if a.OBAtom.GetResidue().GetChain() != lig_obj.chain] if ligand.type == 'INTRA': # Interactions within the chain bs_atoms = [a for a in bs_atoms if a.OBAtom.GetResidue().GetChain() == lig_obj.chain] bs_atoms_refined = [] # Create hash with BSRES -> (MINDIST_TO_LIG, AA_TYPE) # and refine binding site atom selection with exact threshold min_dist = {} for r in bs_atoms: bs_res_id = ''.join([str(whichresnumber(r)), whichchain(r)]) for l in ligand.mol.atoms: distance = euclidean3d(r.coords, l.coords) if bs_res_id not in min_dist: min_dist[bs_res_id] = (distance, whichrestype(r)) elif min_dist[bs_res_id][0] > distance: min_dist[bs_res_id] = (distance, whichrestype(r)) if distance <= config.BS_DIST and r not in bs_atoms_refined: bs_atoms_refined.append(r) num_bs_atoms = len(bs_atoms_refined) logger.info(f'binding site atoms in vicinity ({config.BS_DIST} A max. dist: {num_bs_atoms})') bs_obj = BindingSite(bs_atoms_refined, self.protcomplex, self, self.altconf, min_dist, self.Mapper) pli_obj = PLInteraction(lig_obj, bs_obj, self) self.interaction_sets[ligand.mol.title] = pli_obj def extract_bs(self, cutoff, ligcentroid, resis): """Return list of ids from residues belonging to the binding site""" return [obres.GetIdx() for obres in resis if self.res_belongs_to_bs(obres, cutoff, ligcentroid)] def res_belongs_to_bs(self, res, cutoff, ligcentroid): """Check for each residue if its centroid is within a certain distance to the ligand centroid. Additionally checks if a residue belongs to a chain restricted by the user (e.g. by defining a peptide chain)""" 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) def get_atom(self, idx): return self.atoms[idx] @property def output_path(self): return self._output_path @output_path.setter def output_path(self, path): self._output_path = tilde_expansion(path)