diff --git a/.gitignore b/.gitignore index a193256..1e4d186 100644 --- a/.gitignore +++ b/.gitignore @@ -48,4 +48,6 @@ target/ # Other *~ .nfs* -tests/ \ No newline at end of file +tests/ +# JetBrains IDE settings +.idea/ diff --git a/CHANGES.txt b/CHANGES.txt index 741042a..51b4c85 100644 --- a/CHANGES.txt +++ b/CHANGES.txt @@ -1,5 +1,15 @@ Changelog --------- +# 3.0.0 +* added new `--regions` flag to enable detection of interactions between protein regions (by @snbolz) +* reworked XML and TXT report file naming; changes default behavior +* introduced pyproject.toml to ensure compliance with Python packaging standards +* updated Dockerfile +* fixed switched coordinates of halogen bond acceptor and donor atoms +* fixed handling of PDB input passed as string +* fixed handling of modified residues in protein and peptide ligands +* several minor bug fixes + # 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 diff --git a/plip/basic/config.py b/plip/basic/config.py index d4b67d8..3d6c337 100644 --- a/plip/basic/config.py +++ b/plip/basic/config.py @@ -1,7 +1,7 @@ -__version__ = '2.4.0' +__version__ = '3.0.0' __maintainer__ = 'PharmAI GmbH (2020-2021) - www.pharm.ai - hello@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" +__citation_information__ = "Schake,P. Bolz,SN. et al. PLIP 2025: introducing protein–protein interactions to the protein–ligand interaction profiler. " \ + "Nucl. Acids Res. (10 May 2025), gkaf361. doi: 10.1093/nar/gkaf361" import logging @@ -28,11 +28,13 @@ RESIDUES = {} KEEPMOD = False DNARECEPTOR = False -OUTPUTFILENAME = "report" # Naming for the TXT and XML report files +OUTPUTFILENAME = None # 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 +CHAINS = None # Define chains for protein-protein interaction detection +REGIONS = None +COMPRESS = False # Compress XML and TXT report files # Configuration file for Protein-Ligand Interaction Profiler (PLIP) diff --git a/plip/basic/remote.py b/plip/basic/remote.py index 0c67f6f..89af974 100644 --- a/plip/basic/remote.py +++ b/plip/basic/remote.py @@ -21,11 +21,12 @@ def __init__(self, mol, site): # General Information self.lig_members = sorted(pli.ligand.members) - self.sourcefile = pcomp.sourcefiles['pdbcomplex'] + self.source_pdb_file_content = pcomp.sourcefiles['pdbstring'] # store pdb file content as string self.corrected_pdb = pcomp.corrected_pdb self.pdbid = mol.pymol_name self.hetid = ligand.hetid self.ligandtype = ligand.type + self.regions = ligand.regions self.chain = ligand.chain if not ligand.chain == "0" else "" # #@todo Fix this self.position = str(ligand.position) self.uid = ":".join([self.hetid, self.chain, self.position]) diff --git a/plip/basic/supplemental.py b/plip/basic/supplemental.py index 6939ae5..2f4b990 100644 --- a/plip/basic/supplemental.py +++ b/plip/basic/supplemental.py @@ -55,15 +55,29 @@ def whichchain(atom): return atom.GetResidue().GetChain() if atom.GetResidue() is not None else None -def residue_belongs_to_receptor(res, config): +def residue_belongs_to_receptor(res, regions=None): """tests whether the residue is defined as receptor and is not part of a peptide or residue ligand.""" - if config.CHAINS: + if regions: + ligand_region, bs_region = regions + chain = res.GetChain() + num = res.GetNum() + if ligand_region and bs_region: # both ligand and receptor region were defined + if num in ligand_region.get(chain, []): + return False # residue belongs to ligand + if num in bs_region.get(chain, []): + return True # residue belongs to receptor and not to ligand + return False + else: + return num not in ligand_region.get(chain, []) + + elif 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. @@ -184,9 +198,17 @@ def cluster_doubles(double_list): # File operations ################# -def tilde_expansion(folder_path): +def tilde_expansion(folder_paths): """Tilde expansion, i.e. converts '~' in paths into .""" - return os.path.expanduser(folder_path) if '~' in folder_path else folder_path + if isinstance(folder_paths, list): + expanded_paths = [] + for p in folder_paths: + if "~" in p: + p = os.path.expanduser(p) + expanded_paths.append(p) + return expanded_paths + else: + return os.path.expanduser(folder_paths) if "~" in folder_paths else folder_paths def folder_exists(folder_path): @@ -231,6 +253,18 @@ def start_pymol(quiet=False, options='-p', run=False): pymol.cmd.feedback('disable', 'all', 'everything') +def select_region(region): + """region is a dictionary with chains as keys and a list of residue numbers as values.""" + selection = f"(chain " + for i, (chain, res_numbers) in enumerate(region.items()): + selection += f"{chain} and not resn HOH and resi {'+'.join(map(str, res_numbers))}" + if not i == len(region) - 1: + selection += ") or (chain " + else: + selection += ")" + return selection + + def nucleotide_linkage(residues): """Support for DNA/RNA ligands by finding missing covalent linkages to stitch DNA/RNA together.""" diff --git a/plip/exchange/report.py b/plip/exchange/report.py index 8c489c8..bfb89e4 100644 --- a/plip/exchange/report.py +++ b/plip/exchange/report.py @@ -2,6 +2,7 @@ from operator import itemgetter import lxml.etree as et +import gzip from plip.basic import config from plip.basic.config import __version__ @@ -91,20 +92,31 @@ def get_bindingsite_data(self): else: self.txtreport.append('No interactions detected.') - def write_xml(self, as_string=False): + def write_xml(self, as_string: bool = False): """Write the XML report""" if not as_string: - et.ElementTree(self.xmlreport).write('{}/{}.xml'.format(self.outpath, self.outputprefix), pretty_print=True, - xml_declaration=True) + tree = et.ElementTree(self.xmlreport) + if config.COMPRESS: + with gzip.open(f"{self.outpath}{self.outputprefix}.xml.gz", "wb") as xml_file: + tree.write(xml_file, pretty_print=True, xml_declaration=True, encoding="utf-8") + else: + tree.write(f"{self.outpath}{self.outputprefix}.xml", pretty_print=True, xml_declaration=True, + encoding="utf-8") else: output = et.tostring(self.xmlreport, pretty_print=True) print(output.decode('utf8')) - def write_txt(self, as_string=False): + def write_txt(self, as_string: bool = False): """Write the TXT report""" if not as_string: - with open('{}/{}.txt'.format(self.outpath, self.outputprefix), 'w') as f: - [f.write(textline + '\n') for textline in self.txtreport] + if config.COMPRESS: + with gzip.open(f"{self.outpath}{self.outputprefix}.txt.gz", "wb") as txt_file: + for textline in self.txtreport: + txt_file.write((textline + '\n').encode('utf-8')) + else: + with open(f"{self.outpath}{self.outputprefix}.txt", 'w') as txt_file: + for textline in self.txtreport: + txt_file.write(textline + '\n') else: output = '\n'.join(self.txtreport) print(output) @@ -280,7 +292,7 @@ def __init__(self, plcomplex): '%.2f' % halogen.distance, '%.2f' % halogen.don_angle, '%.2f' % halogen.acc_angle, halogen.don_orig_idx, halogen.donortype, halogen.acc_orig_idx, halogen.acctype, - halogen.acc.o.coords, halogen.don.x.coords)) + halogen.don.x.coords, halogen.acc.o.coords)) ################### # METAL COMPLEXES # diff --git a/plip/plipcmd.py b/plip/plipcmd.py index b79b7df..6f44215 100644 --- a/plip/plipcmd.py +++ b/plip/plipcmd.py @@ -39,18 +39,28 @@ 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'): +def parse_report_filename(parser, name_config): + if name_config is not None: + dir_part, name_config = os.path.split(name_config) + if not name_config: # provided filename is a directory. + parser.error(f"Report filename must be a file basename not a directory.") + if dir_part: # provided filename contains a directory, but file will be written to outpath. + logger.warning(f"Report will be written to {config.OUTPATH}") + base_extensions = name_config.split(".") + name_config = base_extensions[0] # remove all file extensions + + if len(base_extensions) > 1: # Print a warning when improper file extensions were given. + first_ext_invalid = base_extensions[1] not in ["xml", "txt"] + non_compressed_warning = not config.COMPRESS and (len(base_extensions) > 2 or first_ext_invalid) + compressed_warning = len(base_extensions) == 2 or (len(base_extensions) >= 3 and ( + len(base_extensions) > 3 or (first_ext_invalid or base_extensions[2] != "gz"))) + if non_compressed_warning or compressed_warning: + logger.warning("Improper report filename extension(s) will be replaced.") + return name_config + + +def process_pdb(pdbfile, outpath, as_string=False, batch_idx=None): """Analysis of a single PDB file with optional chain filtering.""" if not as_string: pdb_file_name = pdbfile.split('/')[-1] @@ -66,9 +76,6 @@ def process_pdb(pdbfile, outpath, as_string=False, outputprefix='report'): create_folder_if_not_exists(outpath) - # Generate the report files - streport = StructureReport(mol, outputprefix=outputprefix) - config.MAXTHREADS = min(config.MAXTHREADS, len(mol.interaction_sets)) ###################################### @@ -86,11 +93,20 @@ def process_pdb(pdbfile, outpath, as_string=False, outputprefix='report'): else: [visualize_in_pymol(plcomplex) for plcomplex in complexes] - if config.XML: # Generate report in xml format - streport.write_xml(as_string=config.STDOUT) + # Generate the report files + if config.XML or config.TXT: + # set default filename prefix + name = mol.pymol_name.upper() if as_string else os.path.splitext(os.path.basename(pdbfile))[0] + # for batch processing add batch index to custom report names + idx_str = f"_{batch_idx}" if batch_idx is not None else "" + outprefix = f"{name}_report" if config.OUTPUTFILENAME is None else f"{config.OUTPUTFILENAME}{idx_str}" + streport = StructureReport(mol, outputprefix=outprefix) - if config.TXT: # Generate report in txt (rst) format - streport.write_txt(as_string=config.STDOUT) + if config.XML: # Generate report in xml format + streport.write_xml(as_string=config.STDOUT) + + if config.TXT: # Generate report in txt (rst) format + streport.write_txt(as_string=config.STDOUT) def download_structure(inputpdbid): @@ -126,21 +142,21 @@ def remove_duplicates(slist): return unique -def run_analysis(inputstructs, inputpdbids, chains=None): +def run_analysis(inputstructs, inputpdbids): """Main function. Calls functions for processing, report generation and visualization.""" pdbid, pdbpath = None, None + batch_idx = None # @todo For multiprocessing, implement better stacktracing for errors # Print title and version logger.info(f'Protein-Ligand Interaction Profiler (PLIP) {__version__}') logger.info(f'brought to you by: {config.__maintainer__}') logger.info(f'please cite: {config.__citation_information__}') - output_prefix = config.OUTPUTFILENAME if inputstructs is not None: # Process PDB file(s) - num_structures = len(inputstructs) # @question: how can it become more than one file? The tilde_expansion function does not consider this case. + num_structures = len(inputstructs) inputstructs = remove_duplicates(inputstructs) read_from_stdin = False - for inputstruct in inputstructs: + for idx, inputstruct in enumerate(inputstructs): if inputstruct == '-': # @expl: when user gives '-' as input, pdb file is read from stdin inputstruct = sys.stdin.read() read_from_stdin = True @@ -154,19 +170,16 @@ def run_analysis(inputstructs, inputpdbids, chains=None): logger.error('empty PDB file') sys.exit(1) if num_structures > 1: - basename = inputstruct.split('.')[-2].split('/')[-1] - config.OUTPATH = '/'.join([config.BASEPATH, basename]) - output_prefix = 'report' - process_pdb(inputstruct, config.OUTPATH, as_string=read_from_stdin, outputprefix=output_prefix) + batch_idx = idx + process_pdb(inputstruct, config.OUTPATH, as_string=read_from_stdin, batch_idx=batch_idx) else: # Try to fetch the current PDB structure(s) directly from the RCBS server num_pdbids = len(inputpdbids) inputpdbids = remove_duplicates(inputpdbids) - for inputpdbid in inputpdbids: + for idx, inputpdbid in enumerate(inputpdbids): pdbpath, pdbid = download_structure(inputpdbid) if num_pdbids > 1: - config.OUTPATH = '/'.join([config.BASEPATH, pdbid[1:3].upper(), pdbid.upper()]) - output_prefix = 'report' - process_pdb(pdbpath, config.OUTPATH, outputprefix=output_prefix) + batch_idx = idx + process_pdb(pdbpath, config.OUTPATH, batch_idx=batch_idx) if (pdbid is not None or inputstructs is not None) and config.BASEPATH is not None: if config.BASEPATH in ['.', './']: @@ -197,6 +210,10 @@ def main(): action="store_true") parser.add_argument("-t", "--txt", dest="txt", default=False, help="Generate report file in TXT (RST) format", action="store_true") + parser.add_argument("-z", "--gzip", dest="compress", default=False, + help="XML and TXT report files will be gzip compressed.", action="store_true") + parser.add_argument("--name", dest="outputfilename", default=None, + help="Set a filename for the report TXT and XML files.") parser.add_argument("-y", "--pymol", dest="pymol", default=False, help="Additional PyMOL session files", action="store_true") parser.add_argument("--maxthreads", dest="maxthreads", default=multiprocessing.cpu_count(), @@ -221,16 +238,11 @@ def main(): parser.add_argument("--dnareceptor", dest="dnareceptor", default=False, help="Treat nucleic acids as part of the receptor structure (together with any present protein) instead of as a ligand.", action="store_true") - parser.add_argument("--name", dest="outputfilename", default="report", - help="Set a filename for the report TXT and XML files. Will only work when processing single structures.") ligandtype = parser.add_mutually_exclusive_group() # Either peptide/inter or intra mode ligandtype.add_argument("--peptides", "--inter", dest="peptides", default=[], 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") @@ -256,17 +268,19 @@ def main(): help=argparse.SUPPRESS) # Add argument to define receptor and ligand chains - parser.add_argument("--chains", dest="chains", type=str, + ligandtype.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.") + # Add argument to define receptor and ligand regions of the protein + ligandtype.add_argument("--regions", dest="regions", type=str, + help="Specify protein regions as receptor/ligand groups by chain and residue numbers. " + "e.g., ({A: 1-20, 25, 26, 28}, {A: 54-63, B: 17-46}) to define residues" + "1-20, 25, 26, and 28 of chain A as ligand and residues 54-63 of chain A and residues" + "17-46 of chain B as receptor. Defining a receptor is optional. Multiple ligand-receptor" + "pairs can be parsed as a list of tuples [({A: 1-20}, {B: 17-46}), ({C: 5-43}, {D: 7})].") 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 @@ -282,6 +296,7 @@ def main(): config.MAXTHREADS = arguments.maxthreads config.XML = arguments.xml config.TXT = arguments.txt + config.COMPRESS = arguments.compress config.PICS = arguments.pics config.PYMOL = arguments.pymol config.STDOUT = arguments.stdout @@ -293,7 +308,6 @@ 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 @@ -301,9 +315,62 @@ def main(): config.KEEPMOD = arguments.keepmod config.DNARECEPTOR = arguments.dnareceptor config.OUTPUTFILENAME = arguments.outputfilename + if config.OUTPUTFILENAME is not None: + if config.XML or config.TXT: + config.OUTPUTFILENAME = parse_report_filename(parser, config.OUTPUTFILENAME) config.NOHYDRO = arguments.nohydro config.MODEL = arguments.model + def expand_ranges(residue_ranges): + """ + Takes '1-3, 5, 7' -> [1, 2, 3, 5, 7] + """ + parts = [p.strip() for p in residue_ranges.split(',')] + res_list = [] + for p in parts: + if '-' in p: + start, end = p.split('-') + res_list.extend(range(int(start), int(end) + 1)) + else: + res_list.append(int(p)) + return res_list + + try: + # add inner quotes for Python backend and expand residue ranges to lists of residue numbers. + if not arguments.regions: + config.REGIONS = None + else: + import re + # add quotes around keys (chain IDs) + quoted = re.sub(pattern=r'([{,]\s*)([a-zA-Z0-9_]+)\s*:', repl=r'\1"\2":', string=arguments.regions) + # add quotes around values (residue numbers) + quoted = re.sub(pattern=r':\s*([0-9,\s\-]+?)\s*(?=,\s*"[a-zA-Z0-9_]+":|})', repl=r': "\1"', string=quoted) + # add comma to tuples with only one dictionary (without comma would not be treated as tuple) + ensure_tuples = re.sub(r'\((\s*{[^{}]*}\s*)\)', r'(\1,)', quoted) + config.REGIONS = ast.literal_eval(ensure_tuples) # convert string to tuple(s) of one or two dictionaries + if not isinstance(config.REGIONS, list): + # embed single tuple in a list to give the regions config a consistent data structure + config.REGIONS = [config.REGIONS] + if not all(isinstance(lig_rec, tuple) for lig_rec in config.REGIONS) or any( + len(lig_rec) > 2 for lig_rec in config.REGIONS) or not all( + isinstance(reg, dict) for lig_rec in config.REGIONS for reg in lig_rec): + raise ValueError( + "Regions must be specified by tuples of one or two dictionaries (ligand, receptor)") + config.REGIONS = [ + tuple( + {chain: expand_ranges(residues) for chain, residues in dct.items()} + for dct in lig_rec + ) + for lig_rec in config.REGIONS + ] + config.REGIONS = [ + (lig_rec[0], None) if len(lig_rec) == 1 else (lig_rec[0], lig_rec[1]) + for lig_rec in config.REGIONS + ] + except (ValueError, SyntaxError): + parser.error("The --regions option must be in the format '({A: 1-20, 25, 27}, {A: 54-63, B: 17-46})' Multiple region tuples can be given in a list.") + + try: # add inner quotes for python backend if not arguments.chains: @@ -312,7 +379,6 @@ def main(): import re quoted_input = re.sub(r'(? 0: corrected_pdb += ''.join(model_dict[config.MODEL]) corrected_lines += model_dict[config.MODEL] except KeyError: corrected_pdb = ''.join(model_dict[1]) + self.pdb_file_was_corrected = True corrected_lines = model_dict[1] config.MODEL = 1 logger.warning('invalid model number specified, using first model instead') @@ -146,7 +148,10 @@ def fix_pdbline(self, pdbline, lastnum): return None, lastnum # TER Entries also have continuing numbering, consider them as well if pdbline.startswith('TER'): - new_num = lastnum + 1 + if not pdbline[6:11]: # pdb files saved from PyMol skip the number in TER entries + new_num = lastnum + else: + new_num = lastnum + 1 if pdbline.startswith('ATOM'): new_num = lastnum + 1 current_num = int(pdbline[6:11]) @@ -243,8 +248,14 @@ def getpeptides(self, chain): 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 and not self.is_het_residue(o)] # All residues from chain + # All residues from chain + if config.KEEPMOD: + all_from_chain = [o for o in pybel.ob.OBResidueIter( + self.proteincomplex.OBMol) if o.GetChain() == chain and (not self.is_het_residue(o) or + o.GetName() in self.modresidues)] + else: + all_from_chain = [o for o in pybel.ob.OBResidueIter( + self.proteincomplex.OBMol) if o.GetChain() == chain and not self.is_het_residue(o)] if len(all_from_chain) == 0: return None else: @@ -252,12 +263,41 @@ def getpeptides(self, chain): ligand = self.extract_ligand(non_water) return ligand + def getregion(self, ligand_region, bs_region=None): + all_from_region = [] + for chain, residue_numbers in ligand_region.items(): + if config.KEEPMOD: + residues = [ + res for res in pybel.ob.OBResidueIter( + self.proteincomplex.OBMol) if ( + res.GetChain() == chain and + res.GetNum() in residue_numbers and + (not self.is_het_residue(res) or res.GetName() in self.modresidues) + ) + ] + else: + residues = [ + res for res in pybel.ob.OBResidueIter( + self.proteincomplex.OBMol) if ( + res.GetChain() == chain and + res.GetNum() in residue_numbers and not + self.is_het_residue(res) + ) + ] + all_from_region.extend(residues) + if len(all_from_region) == 0: + return None + else: + non_water = [res for res in all_from_region if not res.GetResidueProperty(9)] # 9 is water + ligand = self.extract_ligand(non_water, regions=(ligand_region, bs_region)) + 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 and config.CHAINS is None: + if config.PEPTIDES == [] and config.INTRA is None and config.CHAINS is None and config.REGIONS is None: # Extract small molecule ligands (default) ligands = [] @@ -285,11 +325,17 @@ def getligs(self): 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 and not config.CHAINS: + if config.PEPTIDES: 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.REGIONS: + # regions are defined by tuples of dictionaries in which first dictionary defines the ligand + # and second dictionary defines the receptor + peptide_ligands = [] + for lig_rec in config.REGIONS: + peptide_ligands.append(self.getregion(lig_rec[0], lig_rec[1])) 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 @@ -303,9 +349,9 @@ def getligs(self): return [lig for lig in ligands if len(lig.mol.atoms) != 0] - def extract_ligand(self, kmer): + def extract_ligand(self, kmer, regions=None): """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') + data = namedtuple('ligand', 'mol hetid chain position water members longname type atomorder can_to_pdb regions') 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] @@ -313,7 +359,7 @@ def extract_ligand(self, kmer): names = [x[0] for x in members] longname = '-'.join([x[0] for x in members]) - if config.PEPTIDES or config.CHAINS: + if config.PEPTIDES or config.CHAINS or config.REGIONS: ligtype = 'PEPTIDE' elif config.INTRA is not None: ligtype = 'INTRA' @@ -363,6 +409,10 @@ def extract_ligand(self, kmer): lig.title = ':'.join((rname, rchain, str(rnum))) self.mapper.ligandmaps[lig.title] = mapold + if config.REGIONS: + lig_idx = config.REGIONS.index(regions) + lig.title = ':'.join((rname, rchain, str(rnum), str(lig_idx))) + logger.debug('renumerated molecule generated') if not config.NOPDBCANMAP: @@ -376,7 +426,7 @@ def extract_ligand(self, kmer): 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) + can_to_pdb=can_to_pdb, regions=regions) return ligand @staticmethod @@ -585,6 +635,94 @@ def find_rings(self, mol, all_atoms): type=ring_type)) return rings + def append_func_group_to_data(self, a, a_set, data, a_orig_idx, charge_type, center, fgroup, res=None): + """Appends atoms that are part of a functional group as named tuple to the a_set""" + if not res: + if not isinstance(a, list): + a_orig = self.Mapper.id_to_atom(a_orig_idx) + a = [a, ] + a_orig = [a_orig, ] + a_orig_idx = [a_orig_idx, ] + else: + a_orig = [self.Mapper.id_to_atom(idx) for idx in a_orig_idx] + a_set.append(data(atoms=a, orig_atoms=a_orig, atoms_orig_idx=a_orig_idx, type=charge_type, + center=center, fgroup=fgroup)) + return a_set + else: + if not isinstance(a, list): + a = [a, ] + a_orig_idx = [a_orig_idx, ] + a_set.append(data(atoms=a, atoms_orig_idx=a_orig_idx, type=charge_type, center=center, + restype=res.GetName(), resnr=res.GetNum(), reschain=res.GetChain())) + return a_set + + def append_if_charged_func_group(self, a, a_set, data, res=None): + """Checks if atom is part of a charged functional group and appends it to a_set if True.""" + a_orig_idx = self.Mapper.mapid(a.idx, mtype=self.mtype, bsid=self.bsid) + if self.is_functional_group(a, 'quartamine'): + a_set = self.append_func_group_to_data(a=a, a_set=a_set, data=data, a_orig_idx=a_orig_idx, + charge_type='positive', center=list(a.coords), fgroup='quartamine', + res=res) + elif self.is_functional_group(a, 'tertamine'): + a_set = self.append_func_group_to_data(a=a, a_set=a_set, data=data, a_orig_idx=a_orig_idx, + charge_type='positive', center=list(a.coords), fgroup='tertamine', + res=res) + if self.is_functional_group(a, 'sulfonium'): + a_set = self.append_func_group_to_data(a=a, a_set=a_set, data=data, a_orig_idx=a_orig_idx, + charge_type='positive', center=list(a.coords), fgroup='sulfonium', + res=res) + 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] + a_set = self.append_func_group_to_data(a=a_contributing, a_set=a_set, data=data, + a_orig_idx=a_contributing_orig_idx, + charge_type='negative', center=a.coords, fgroup='phosphate', + res=res) + 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] + a_set = self.append_func_group_to_data(a=a_contributing, a_set=a_set, data=data, + a_orig_idx=a_contributing_orig_idx, + charge_type='negative', center=a.coords, fgroup='sulfonicacid', + res=res) + 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)] + a_set = self.append_func_group_to_data(a=a_contributing, a_set=a_set, data=data, + a_orig_idx=a_contributing_orig_idx, + charge_type='negative', center=a.coords, fgroup='sulfate', + res=res) + 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] + a_set = self.append_func_group_to_data(a=a_contributing, a_set=a_set, data=data, + a_orig_idx=a_contributing_orig_idx, + charge_type='negative', + center=centroid([a.coords for a in a_contributing]), + fgroup='carboxylate', res=res) + 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] + a_set = self.append_func_group_to_data(a=a_contributing, a_set=a_set, data=data, + a_orig_idx=a_contributing_orig_idx, + charge_type='positive', center=a.coords, fgroup='guanidine', + res=res) + return a_set + def get_hydrophobic_atoms(self): return self.hydroph_atoms @@ -603,6 +741,49 @@ def get_pos_charged(self): def get_neg_charged(self): return [charge for charge in self.charged if charge.type == 'negative'] + @staticmethod + def is_functional_group(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 + class PLInteraction: """Class to store a ligand, a protein and their interactions.""" @@ -752,7 +933,7 @@ def refine_hydrophobic(all_h, pistacks): hydroph = [h for h in sel2.values()] hydroph_final = [] # 3. If a protein atom interacts with several neighboring ligand atoms, just keep the one with the closest dist - if config.PEPTIDES or config.INTRA or config.CHAINS: + if config.PEPTIDES or config.INTRA or config.CHAINS or config.REGIONS: # the ligand also consists of amino acid residues, repeat step 2 just the other way around sel3 = {} for h in hydroph: @@ -925,13 +1106,14 @@ def refine_water_bridges(wbridges, hbonds_ldon, hbonds_pdon): class BindingSite(Mol): - def __init__(self, atoms, protcomplex, cclass, altconf, min_dist, mapper): + def __init__(self, atoms, protcomplex, cclass, altconf, min_dist, mapper, regions): """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.regions = regions 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) @@ -960,7 +1142,7 @@ def find_charged(self, 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 residue_belongs_to_receptor(r, config)]: + for res in [r for r in pybel.ob.OBResidueIter(mol.OBMol) if residue_belongs_to_receptor(r, self.regions)]: if config.INTRA is not None: if res.GetChain() != config.INTRA: continue @@ -980,7 +1162,7 @@ def find_charged(self, mol): restype=res.GetName(), resnr=res.GetNum(), reschain=res.GetChain())) - if res.GetName() in ('GLU', 'ASP'): # Aspartic or Glutamic Acid + elif 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: @@ -1001,10 +1183,15 @@ def find_charged(self, 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())) + if config.KEEPMOD and res.GetName() in self.complex.modres: + atom_indices = [a.GetIdx() for a in pybel.ob.OBResidueAtomIter(res)] + atoms = [atm for atm in self.all_atoms if atm.idx in atom_indices] + for a in atoms: + a_set = self.append_if_charged_func_group(a=a, a_set=a_set, data=data, res=res) return a_set def find_metal_binding(self, mol): @@ -1060,12 +1247,12 @@ def __init__(self, cclass, ligand): self.complex = cclass self.molecule = ligand.mol # Pybel Molecule # 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.smiles = "" if (config.INTRA or config.PEPTIDES or config.CHAINS or config.REGIONS) 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: self.smiles = self.smiles.split()[0] - else: + elif not (config.INTRA or config.PEPTIDES or config.CHAINS or config.REGIONS): logger.warning(f'could not write SMILES for ligand {ligand}') self.smiles = '' self.heavy_atoms = self.molecule.OBMol.NumHvyAtoms() # Heavy atoms count @@ -1081,6 +1268,7 @@ def __init__(self, cclass, ligand): self.molweight, self.logp = float(descvalues['MW']), float(descvalues['logP']) self.num_rot_bonds = int(self.molecule.OBMol.NumRotors()) self.atomorder = ligand.atomorder + self.regions = ligand.regions ########################################################## # Special Case for hydrogen bond acceptor identification # @@ -1134,49 +1322,6 @@ 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] - @staticmethod - def is_functional_group(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') @@ -1201,75 +1346,9 @@ def find_charged(self, all_atoms): """ data = namedtuple('lcharge', 'atoms orig_atoms atoms_orig_idx type center fgroup') a_set = [] - if not (config.INTRA or config.PEPTIDES or config.CHAINS): + if not (config.INTRA or config.PEPTIDES or config.CHAINS or config.REGIONS): 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')) + a_set = self.append_if_charged_func_group(a=a, a_set=a_set, data=data) else: """We have peptide/protein chain as ligand""" """Looks for positive charges in arginine, histidine or lysine, for negative in aspartic and glutamic acid.""" @@ -1292,7 +1371,7 @@ def find_charged(self, all_atoms): 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 + elif 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: @@ -1305,6 +1384,11 @@ def find_charged(self, all_atoms): type='negative', center=centroid([ac.coords for ac in a_contributing]), fgroup=res.GetName()+str(res.GetNum())+res.GetChain())) + if config.KEEPMOD and res.GetName() in self.complex.modres: + atom_indices = [a.GetIdx() for a in pybel.ob.OBResidueAtomIter(res)] + atoms = [atm for atm in all_atoms if atm.idx in atom_indices] + for a in atoms: + a_set = self.append_if_charged_func_group(a=a, a_set=a_set, data=data) return a_set def find_metal_binding(self, lig_atoms, water_oxygens): @@ -1446,7 +1530,9 @@ def load_pdb(self, pdbpath, as_string=False): if not as_string: self.sourcefiles['filename'] = os.path.basename(self.sourcefiles['pdbcomplex']) - self.protcomplex, self.filetype = read_pdb(self.corrected_pdb, as_string= (self.corrected_pdb != pdbpath)) # self.corrected_pdb may fallback to pdbpath + self.protcomplex, self.filetype = read_pdb(self.corrected_pdb, as_string= (self.corrected_pdb != pdbpath)) # self.corrected_pdb may fallback to pdbpath + else: + 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 @@ -1487,11 +1573,19 @@ def load_pdb(self, pdbpath, as_string=False): if len(self.excluded) != 0: logger.info(f'excluded molecules as ligands: {self.excluded}') - if config.DNARECEPTOR: + if config.DNARECEPTOR and config.KEEPMOD: self.resis = [obres for obres in pybel.ob.OBResidueIter( self.protcomplex.OBMol) if obres.GetName() in config.DNA + config.RNA ] + [obres for obres in pybel.ob.OBResidueIter( - self.protcomplex.OBMol) if obres.GetResidueProperty(0)] + self.protcomplex.OBMol) if obres.GetResidueProperty(0) or obres.GetName() in self.modres] + elif config.DNARECEPTOR: + self.resis = [obres for obres in pybel.ob.OBResidueIter( + self.protcomplex.OBMol) if obres.GetName() in config.DNA + config.RNA + ] + [obres for obres in pybel.ob.OBResidueIter( + self.protcomplex.OBMol) if obres.GetResidueProperty(0)] + elif config.KEEPMOD: + self.resis = [obres for obres in pybel.ob.OBResidueIter( + self.protcomplex.OBMol) if obres.GetResidueProperty(0) or obres.GetName() in self.modres] else: self.resis = [obres for obres in pybel.ob.OBResidueIter( self.protcomplex.OBMol) if obres.GetResidueProperty(0)] @@ -1504,6 +1598,11 @@ def load_pdb(self, pdbpath, as_string=False): else: logger.info(f'structure contains no ligands') + if as_string or pdbparser.pdb_file_was_corrected: + self.sourcefiles['pdbstring'] = self.corrected_pdb + else: + self.sourcefiles['pdbstring'] = open(pdbpath, 'r').read() + def analyze(self): """Triggers analysis of all complexes in structure""" for ligand in self.ligands: @@ -1532,12 +1631,17 @@ def characterize_complex(self, ligand): 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) + + resis = self.resis + if config.KEEPMOD: + resis = self.exclude_ligand_modresidues(lig_obj.members, resis) + + bs_res = self.extract_bs(cutoff, lig_obj.centroid, resis, lig_obj.regions) # 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 ligand.type == 'PEPTIDE' and not config.REGIONS: # 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': @@ -1561,16 +1665,25 @@ def characterize_complex(self, ligand): 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) + bs_obj = BindingSite(bs_atoms_refined, self.protcomplex, self, self.altconf, min_dist, self.Mapper, lig_obj.regions) pli_obj = PLInteraction(lig_obj, bs_obj, self) self.interaction_sets[ligand.mol.title] = pli_obj - def extract_bs(self, cutoff, ligcentroid, resis): + def exclude_ligand_modresidues(self, ligmembers, resis): + """If the ligand contains modified residues, exclude these from the receptor residues.""" + lig_modres = [member for member in ligmembers if member[0] in self.modres] + if lig_modres: + return [obres for obres in resis if + not (obres.GetName(), obres.GetChain(), int32_to_negative(obres.GetNum())) in lig_modres] + else: + return resis + + def extract_bs(self, cutoff, ligcentroid, resis, regions=None): """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)] + return [obres.GetIdx() for obres in resis if self.res_belongs_to_bs(obres, cutoff, ligcentroid, regions)] @staticmethod - def res_belongs_to_bs(res, cutoff, ligcentroid): + def res_belongs_to_bs(res, cutoff, ligcentroid, regions=None): """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)]) @@ -1578,7 +1691,7 @@ def res_belongs_to_bs(res, cutoff, ligcentroid): near_enough = True if euclidean3d(rescentroid, ligcentroid) < cutoff else False #Todo: Test if properly working # Add restriction via chains flag - return near_enough and residue_belongs_to_receptor(res, config) + return near_enough and residue_belongs_to_receptor(res, regions) def get_atom(self, idx): return self.atoms[idx] diff --git a/plip/visualization/visualize.py b/plip/visualization/visualize.py index b419500..19fd71b 100644 --- a/plip/visualization/visualize.py +++ b/plip/visualization/visualize.py @@ -1,7 +1,7 @@ from pymol import cmd from plip.basic import config, logger -from plip.basic.supplemental import start_pymol +from plip.basic.supplemental import start_pymol, select_region from plip.visualization.pymol import PyMOLVisualizer logger = logger.get_logger() @@ -21,6 +21,9 @@ def visualize_in_pymol(plcomplex): chain = plcomplex.chain if config.PEPTIDES or config.CHAINS: vis.ligname = 'PeptideChain%s' % plcomplex.chain + if config.REGIONS: + lig_idx = config.REGIONS.index(plcomplex.regions) + vis.ligname = f'Ligand{lig_idx+1}' if config.INTRA is not None: vis.ligname = 'Intra%s' % plcomplex.chain @@ -37,7 +40,8 @@ def visualize_in_pymol(plcomplex): start_pymol(run=True, options='-pcq', quiet=not config.VERBOSE and not config.SILENT) vis.set_initial_representations() - cmd.load(plcomplex.sourcefile) + cmd.read_pdbstr(plcomplex.source_pdb_file_content, f"{pdbid}_orig") + cmd.frame(config.MODEL) current_name = cmd.get_object_list(selection='(all)')[0] @@ -45,10 +49,11 @@ def visualize_in_pymol(plcomplex): cmd.set_name(current_name, pdbid) cmd.hide('everything', 'all') if config.PEPTIDES: - 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) + cmd.select(ligname, 'chain %s and not resn HOH' % plcomplex.chain) + elif config.REGIONS: + lig_region, rec_region = plcomplex.regions + lig_selection = select_region(lig_region) + cmd.select(ligname, lig_selection) 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}') @@ -104,6 +109,10 @@ def visualize_in_pymol(plcomplex): filename = "%s_PeptideChain%s" % (pdbid.upper(), plcomplex.chain) if config.PYMOL: vis.save_session(config.OUTPATH, override=filename) + elif config.REGIONS: + filename = f"{pdbid.upper()}_ProteinRegion_{vis.ligname}" + if config.PYMOL: + vis.save_session(config.OUTPATH, override=filename) elif config.INTRA is not None: filename = "%s_IntraChain%s" % (pdbid.upper(), plcomplex.chain) if config.PYMOL: diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..a910db7 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,54 @@ +[build-system] +build-backend = "setuptools.build_meta" +requires = ["setuptools","pip","requests"] + +[project] +name = "plip" +license = "GPL-2.0-only" +dynamic = ["version"] +description = "PLIP - Fully automated protein-ligand interaction profiler" +readme = "README.md" +requires-python = ">=3.6" + +authors = [ + {name = "PharmAI GmbH", email = "hello@pharm.ai"} +] +maintainers = [ + {name = "PharmAI GmbH", email = "hello@pharm.ai"} +] + +keywords = ["bioinformatics", "protein-ligand", "interactions", "molecular-modeling"] + +classifiers = [ + "Development Status :: 5 - Production/Stable", + "Intended Audience :: Science/Research", + "Natural Language :: English", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.6", + "Programming Language :: Python :: 3.7", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Topic :: Scientific/Engineering :: Bio-Informatics" +] + +dependencies = [ + "numpy>=1.13.3", + "lxml>=4.2.1", + "openbabel>=3.1.1", +] + +[project.urls] +Homepage = "https://github.com/pharmai/plip" +Repository = "https://github.com/pharmai/plip" + +[project.scripts] +plip = "plip.plipcmd:main" + +[tool.setuptools.packages.find] +include = ["plip*"] + +[tool.setuptools.dynamic] +version = {attr = "plip.basic.config.__version__"} diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index 6d4d07c..0000000 --- a/requirements.txt +++ /dev/null @@ -1,4 +0,0 @@ -lxml~=4.2.1 -numpy~=1.13.3 -pymol~=2.3.0 -openbabel~=3.0.0 \ No newline at end of file diff --git a/setup.py b/setup.py index fb2404b..d5cbdec 100644 --- a/setup.py +++ b/setup.py @@ -3,8 +3,6 @@ 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 @@ -74,7 +72,7 @@ def run(self): return setup(name='plip', - version=config.__version__, + version='3.0.0', description='PLIP - Fully automated protein-ligand interaction profiler', classifiers=[ 'Development Status :: 5 - Production/Stable',