import ost from ost import io from ost import conop from ost import mol from ost import seq def CleanHydrogens(ent, clib): """ Scoring helper - Returns copy of *ent* without hydrogens Non-standard hydrogen naming can cause trouble in residue property assignment which is done by the :class:`ost.conop.RuleBasedProcessor` when loading. In fact, residue property assignment is not done for every residue that has unknown atoms according to the chemical component dictionary. This function therefore re-processes the entity after removing hydrogens. :param ent: Entity to clean :type ent: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView` :param clib: Compound library to perform re-processing after hydrogen removal. :type clib: :class:`ost.conop.CompoundLib` :returns: Cleaned and re-processed ent """ cleaned_ent = mol.CreateEntityFromView(ent.Select( "ele != H and ele != D"), include_exlusive_atoms=False) # process again to set missing residue properties due to non standard # hydrogens processor = conop.RuleBasedProcessor(clib) processor.Process(cleaned_ent) return cleaned_ent def MMCIFPrep(mmcif_path, biounit=None, extract_nonpoly=False, fault_tolerant=False, extract_seqres_mapping=False): """ Scoring helper - Prepares input from mmCIF Only performs gentle cleanup of hydrogen atoms. Further cleanup is delegated to scoring classes. Depending on input flags, the following outputs can be retrieved: * poly_ent (:class:`ost.mol.EntityHandle`): An OpenStructure entity with only polymer chains. This is based on _entity.type extracted from *mmcif_path*. If _entity.type is not defined for every chain, a warning is logged and the returned poly_ent is a selection for peptide and nucleotide residues as defined in the chemical component dictionary. * non_poly_entities (:class:`list` of :class:`ost.mol.EntityHandle`): OpenStructure entities representing all non-polymer (ligand) entities. This is based on _entity.type extracted from *mmcif_path*. If _entity.type is not defined for every chain, an error is raised. * seqres (:class:`ost.seq.SequenceList`): Seqres sequences with entity id as sequence names and the respective canonical seqres as sequence. Set to None and triggers a warning if information is missing in *mmcif_path*. * trg_seqres_mapping (:class:`dict`): Dictionary with chain names in poly_ent as keys and the respective entity ids as values. Set to None and triggers a warning if information is missing in *mmcif_path*. :param mmcif_path: Path to mmCIF file that contains polymer and optionally non-polymer entities :type mmcif_path: :class:`str` :param biounit: If given, construct specified biounit from mmCIF AU :type biounit: :class:`str` :param extract_nonpoly: Controls return value :type extract_nonpoly: :class:`bool` :param fault_tolerant: Passed as parameter to :func:`ost.io.LoadMMCIF` :type fault_tolerant: :class:`bool` :param extract_seqres_mapping: Controls return value :type extract_seqres_mapping: :class:`bool` :returns: poly_ent if *extract_nonpoly*/*extract_seqres_mapping* are False. (poly_ent, non_poly_entities) if *extract_nonpoly* is True. (poly_ent, seqres, trg_seqres_mapping) if *extract_seqres_mapping* is True. (poly_ent, non_poly_entities, seqres, trg_seqres_mapping) if both flags are True. """ clib = conop.GetDefaultLib() if not clib: ost.LogError("A compound library is required. " "Please refer to the OpenStructure website: " "https://openstructure.org/docs/conop/compoundlib/.") raise RuntimeError("No compound library found") # return variables that will be defined depending on input flags poly_ent = None non_poly_entities = None seqres = None trg_seqres_mapping = None # increase loglevel, as we would pollute the info log with weird stuff ost.PushVerbosityLevel(ost.LogLevel.Error) mmcif_entity, mmcif_seqres, mmcif_info = io.LoadMMCIF(mmcif_path, seqres=True, info=True, fault_tolerant=fault_tolerant) # restore old loglevel and return ost.PopVerbosityLevel() mmcif_entity = CleanHydrogens(mmcif_entity, clib) # construct biounit if necessary if biounit is not None: biounit_found = False for bu in mmcif_info.biounits: if bu.id == biounit: mmcif_entity = mol.alg.CreateBU(mmcif_entity, bu) biounit_found = True break if not biounit_found: raise RuntimeError(f"Specified biounit '{biounit}' not in " f"{mmcif_path}") # check if we have entity types defined for each chain missing_entity_types = list() for ch in mmcif_entity.chains: cname = None if biounit is not None: # if a biounit is constructed, you get chain names like: 1.YOLO # we cannot simply split by '.' since '.' is an allowed character # in chain names. => split by first occurence dot_index = ch.name.find('.') if dot_index == -1: cname = ch.name else: cname = ch.name[dot_index+1:] else: cname = ch.name try: entity_id = mmcif_info.GetMMCifEntityIdTr(cname) # the following raises if there is no desc for entity_id entity_desc = mmcif_info.GetEntityDesc(entity_id) except: missing_entity_types.append(cname) if len(missing_entity_types) > 0: msg = f"mmCIF file does not define _entity.type for chains " msg += f"{missing_entity_types}. " if not fault_tolerant: msg += f"Use fault tolerant mode to ignore this error and " msg += f"fallback to select polymers based on peptide/nucleotide " msg += f"residues as defined in chemical component dictionary." raise RuntimeError(msg) msg += f"Fallback to select polymers based on peptide/nucleotide " msg += f"residues as defined in chemical component dictionary (fault " msg += f"tolerant mode)." ost.LogWarning(msg) poly_sel = mmcif_entity.Select("peptide=true or nucleotide=true") poly_ent = mol.CreateEntityFromView(poly_sel, True) else: polymer_entity_ids = mmcif_info.GetEntityIdsOfType("polymer") for ch in mmcif_entity.chains: cname = None if biounit is not None: # if a biounit is constructed, you get chain names like: 1.YOLO # we cannot simply split by '.' since '.' is an allowed character # in chain names. => split by first occurence dot_index = ch.name.find('.') if dot_index == -1: cname = ch.name else: cname = ch.name[dot_index+1:] else: cname = ch.name if mmcif_info.GetMMCifEntityIdTr(cname) in polymer_entity_ids: ch.SetIntProp("poly", 1) poly_sel = mmcif_entity.Select("gcpoly:0=1") poly_ent = mol.CreateEntityFromView(poly_sel, True) if extract_nonpoly: if len(missing_entity_types) > 0: msg = f"mmCIF file does not contain _entity.type for the following " msg += f"chain(s): {missing_entity_types}. Extracting non-polymers " msg += f"from mmCIF files requires _entity.type to be set for all " msg += f"chains." raise RuntimeError(msg) non_polymer_entity_ids = mmcif_info.GetEntityIdsOfType("non-polymer") nonpoly_id = 1 for ch in mmcif_entity.chains: cname = None if biounit is not None: # if a biounit is constructed, you get chain names like: 1.YOLO # we cannot simply split by '.' since '.' is an allowed character # in chain names. => split by first occurence dot_index = ch.name.find('.') if dot_index == -1: cname = ch.name else: cname = ch.name[dot_index+1:] else: cname = ch.name if mmcif_info.GetMMCifEntityIdTr(cname) in non_polymer_entity_ids: ch.SetIntProp("nonpolyid", nonpoly_id) nonpoly_id += 1 non_poly_entities = list() for i in range(1,nonpoly_id): view = mmcif_entity.Select(f"gcnonpolyid:0={i}") if view.GetResidueCount() != 1: raise RuntimeError(f"Expected non-polymer entities in " f"{mmcif_path} to contain exactly 1 " f"residue. Got {view.GetResidueCount()} " f"in chain {view.chains[0].name}") compound = clib.FindCompound(view.residues[0].name) if compound is None: error_msg = f"\"{view.residues[0].name}\" is not available in " \ f"the compound library." if fault_tolerant: error_msg += f"A distance-based heuristic was used to " \ f"connect the ligand atoms (fault tolerant " \ f"mode)." ost.LogWarning(error_msg) else: error_msg += f"Use fault tolerant mode to ignore this " \ f"error and use a distance based heuristic " \ f"to connect the ligand atoms." raise RuntimeError(error_msg) non_poly_entities.append(mol.CreateEntityFromView(view, True)) if extract_seqres_mapping: # mmcif seqres is a list of sequences that relates to # chain names in the assymetric unit. What we want is a list # of sequences that relate to the underlying entities. seqres = seq.CreateSequenceList() seqres_processed = set() for s in mmcif_seqres: entity_id = mmcif_info.GetMMCifEntityIdTr(s.GetName()) if entity_id not in seqres_processed: seqres_processed.add(entity_id) seqres.AddSequence(seq.CreateSequence(entity_id, s.GetGaplessString())) # check if we have SEQRES defined for each polymer chain missing_seqres = list() for ch in poly_ent.chains: cname = None if biounit is not None: # if a biounit is constructed, you get chain names like: 1.YOLO # we cannot simply split by '.' since '.' is an allowed character # in chain names. => split by first occurence dot_index = ch.name.find('.') if dot_index == -1: cname = ch.name else: cname = ch.name[dot_index+1:] else: cname = ch.name entity_id = mmcif_info.GetMMCifEntityIdTr(cname) if entity_id not in seqres_processed: missing_seqres.append(cname) if len(missing_seqres) > 0: msg = f"Extracting chem grouping from mmCIF file requires all " msg += f"SEQRES information set. SEQRES is missing for polymer " msg += f"chain(s) {missing_seqres}. " if not fault_tolerant: msg += f"Use fault tolerant mode to ignore this error and " msg += f"fallback sequence identity-based chem grouping." raise RuntimeError(msg) msg += f"Chem grouping will be based on sequence identity (fault " msg += f"tolerant mode)." ost.LogWarning(msg) seqres = None trg_seqres_mapping = None else: trg_seqres_mapping = dict() if biounit is None: cnames = [ch.name for ch in poly_ent.chains] for cname in cnames: trg_seqres_mapping[cname] = mmcif_info.GetMMCifEntityIdTr(cname) else: bu_cnames = [ch.name for ch in poly_ent.chains] au_cnames = list() for bu_cname in bu_cnames: dot_idx = bu_cname.index(".") au_cnames.append(bu_cname[dot_idx + 1 :]) for au_cname, bu_cname in zip(au_cnames, bu_cnames): trg_seqres_mapping[bu_cname] = mmcif_info.GetMMCifEntityIdTr(au_cname) if extract_nonpoly and extract_seqres_mapping: return (poly_ent, non_poly_entities, seqres, trg_seqres_mapping) elif extract_nonpoly: return (poly_ent, non_poly_entities) elif extract_seqres_mapping: return (poly_ent, seqres, trg_seqres_mapping) else: return poly_ent def PDBPrep(pdb_path, fault_tolerant=False): """ Scoring helper - Prepares scoring input from PDB Only performs gentle cleanup of hydrogen atoms. Further cleanup is delegated to scoring classes. There is no logic to extract ligands from PDB files. Ligands must be provided separately as SDF files in these cases. :param pdb_path: Path to PDB file that contains polymer entities :type pdb_path: :class:`str` :param fault_tolerant: Passed as parameter to :func:`ost.io.LoadPDB` :type fault_tolerant: :class:`bool` :returns: :class:`EntityHandle` from loaded file. """ clib = conop.GetDefaultLib() if not clib: ost.LogError("A compound library is required. " "Please refer to the OpenStructure website: " "https://openstructure.org/docs/conop/compoundlib/.") raise RuntimeError("No compound library found") # increase loglevel, as we would pollute the info log with weird stuff ost.PushVerbosityLevel(ost.LogLevel.Error) pdb_entity = io.LoadPDB(pdb_path, fault_tolerant=fault_tolerant) # restore old loglevel and return ost.PopVerbosityLevel() pdb_entity = CleanHydrogens(pdb_entity, clib) return pdb_entity __all__ = ('CleanHydrogens', 'MMCIFPrep', 'PDBPrep')