From 65ecccca511fc30274cdfd0f02725f05e9c054fd Mon Sep 17 00:00:00 2001 From: Joani Mato Date: Fri, 10 Jul 2026 13:15:49 -0700 Subject: [PATCH] feat: account for radical cofactors --- docs/source/rec_cli_options.rst | 4 ++++ meeko/chemtempgen.py | 33 ++++++++++++++++++++++---------- meeko/cli/mk_prepare_receptor.py | 8 +++++++- meeko/polymer.py | 4 +++- meeko/preparation.py | 3 +++ test/chemtempgen_test.py | 10 ++++++++++ 6 files changed, 50 insertions(+), 12 deletions(-) diff --git a/docs/source/rec_cli_options.rst b/docs/source/rec_cli_options.rst index 0cce471d..3d769c6e 100644 --- a/docs/source/rec_cli_options.rst +++ b/docs/source/rec_cli_options.rst @@ -107,6 +107,10 @@ Receptor Perception Options Default is gasteiger, 'zero' sets all zeros, 'read' requires --read_pqr +.. option:: --radical_cofactor + Meeko fails if the receptor contains a radical cofactor. This option relaxes some of the + checks so that radical cofactors are handled properly. + Grid Box Options ~~~~~~~~~~~~~~~~ diff --git a/meeko/chemtempgen.py b/meeko/chemtempgen.py index 236513f5..b14b812c 100644 --- a/meeko/chemtempgen.py +++ b/meeko/chemtempgen.py @@ -372,7 +372,7 @@ def recharge(rwmol: Chem.RWMol) -> Chem.RWMol: # Attribute Formatters -def get_smiles_with_atom_names(mol: Chem.Mol) -> tuple[str, list[str]]: +def get_smiles_with_atom_names(mol: Chem.Mol, radical_cofactor) -> tuple[str, list[str]]: """ Generate SMILES with atom names in the order of SMILES output @@ -380,6 +380,8 @@ def get_smiles_with_atom_names(mol: Chem.Mol) -> tuple[str, list[str]]: ---------- mol : Chem.Mol Input molecule + radical_cofactor: + true of the input is meant to be a radical Returns ------- @@ -387,7 +389,10 @@ def get_smiles_with_atom_names(mol: Chem.Mol) -> tuple[str, list[str]]: Tuple containing the SMILES string and a list of atom names in the order of SMILES output """ # allHsExplicit may expose the implicit Hs of linker atoms to Smiles; the implicit Hs don't have names - smiles_exh = Chem.MolToSmiles(mol, allHsExplicit=True) + if radical_cofactor: + smiles_exh = Chem.MolToSmiles(mol) + else: + smiles_exh = Chem.MolToSmiles(mol, allHsExplicit=True) smiles_atom_output_order = mol.GetProp('_smilesAtomOutputOrder') delimiters = ('[', ']', ',') @@ -515,7 +520,7 @@ def __eq__(self, other): return False @classmethod - def from_cif(cls, source_cif: str, resname: str): + def from_cif(cls, source_cif: str, resname: str, radical_cofactor): """ Create ChemicalComponent from a chemical component dict file and a searchable residue name in file. @@ -553,6 +558,8 @@ def from_cif(cls, source_cif: str, resname: str): if len(element)==2: element = element[0] + element[1].lower() rdkit_atom = Chem.Atom(element) + if radical_cofactor: + rdkit_atom.SetNoImplicit(True) for attr in atom_attributes: rdkit_atom.SetProp(attr, atom_cols[attr][idx]) # strip double quotes in names @@ -607,6 +614,8 @@ def from_cif(cls, source_cif: str, resname: str): # Finish eidting mol try: rwmol.UpdatePropertyCache() + if radical_cofactor: + Chem.SanitizeMol(rwmol) except Exception as e: logger.error(f"Failed to create rdkitmol from cif. Error: {e} -> template for {resname} will be None. ") return None @@ -620,7 +629,7 @@ def from_cif(cls, source_cif: str, resname: str): rdkit_mol = rwmol.GetMol() # Get Smiles with explicit Hs and ordered atom names - smiles_exh, atom_name = get_smiles_with_atom_names(rdkit_mol) + smiles_exh, atom_name = get_smiles_with_atom_names(rdkit_mol, radical_cofactor) return cls(rdkit_mol, resname, smiles_exh, atom_name) @@ -693,7 +702,7 @@ def make_capped(self, allowed_smarts, capping_names = None, capping_smarts_loc = def make_pretty_smiles(self): """Build and name explicit hydrogens for atoms with implicit Hs by atom names and/or patterns.""" - self.smiles_exh, self.atom_name = get_smiles_with_atom_names(self.rdkit_mol) + self.smiles_exh, self.atom_name = get_smiles_with_atom_names(self.rdkit_mol, False) self.smiles_exh = get_pretty_smiles(self.smiles_exh) return self @@ -912,7 +921,9 @@ def fetch_from_pdb(resname: str, max_retries = 5, backoff_factor = 2, ignore_htt } # Standard pipelines -def build_noncovalent_CC(basename: str, ignore_https_cert = False) -> ChemicalComponent: +def build_noncovalent_CC(basename: str, + ignore_https_cert = False, + radical_cofactor=False) -> ChemicalComponent: """ Build a noncovalent chemical component from a CIF file. @@ -927,7 +938,7 @@ def build_noncovalent_CC(basename: str, ignore_https_cert = False) -> ChemicalCo The constructed ChemicalComponent instance. """ with ChemicalComponent_LoggingControler(): - cc_from_cif = ChemicalComponent.from_cif(fetch_from_pdb(basename, ignore_https_cert=ignore_https_cert), basename) + cc_from_cif = ChemicalComponent.from_cif(fetch_from_pdb(basename, ignore_https_cert=ignore_https_cert), basename, radical_cofactor) if cc_from_cif is None: return None @@ -939,7 +950,8 @@ def build_noncovalent_CC(basename: str, ignore_https_cert = False) -> ChemicalCo err = f"Template Generation failed for {cc.resname}. Error: Molecule breaks into fragments during the deleterious editing. " raise RuntimeError(err) - cc = cc.make_pretty_smiles() + if not radical_cofactor: + cc = cc.make_pretty_smiles() # Check try: @@ -1057,7 +1069,8 @@ class NA_recipe: def build_linked_CCs(basename: str, embed_allowed_smarts: str = None, cap_allowed_smarts: str = None, cap_protonate: bool = False, pattern_to_label_mapping_standard = dict[str, str], - variant_dict = dict[str, tuple], ignore_https_cert = False) -> list[ChemicalComponent]: + variant_dict = dict[str, tuple], ignore_https_cert = False, + radical_cofactor=False) -> list[ChemicalComponent]: """ Build a linked chemical component from a CIF file. @@ -1084,7 +1097,7 @@ def build_linked_CCs(basename: str, embed_allowed_smarts: str = None, List of ChemicalComponent instances with the added variants. """ with ChemicalComponent_LoggingControler(): - cc_from_cif = ChemicalComponent.from_cif(fetch_from_pdb(basename, ignore_https_cert=ignore_https_cert), basename) + cc_from_cif = ChemicalComponent.from_cif(fetch_from_pdb(basename, ignore_https_cert=ignore_https_cert), basename, radical_cofactor) if cc_from_cif is None: return None diff --git a/meeko/cli/mk_prepare_receptor.py b/meeko/cli/mk_prepare_receptor.py index a3373a81..36385a9b 100755 --- a/meeko/cli/mk_prepare_receptor.py +++ b/meeko/cli/mk_prepare_receptor.py @@ -255,7 +255,6 @@ def get_args(): help='specify the residues for which to make terminal functional group rotatable by the chain ID and residue number, e.g. -t ":42,B:23" is equivalent to -t ":42" -t "B:23" (leave chain ID empty if omitted in input PDB or mmCIF)', ) - config_group.add_argument( "--compute_charges", help="compute charges from scratch with the given charge model instead of reading from template (note: this option is slower)", @@ -268,6 +267,11 @@ def get_args(): help="default is gasteiger, 'zero' sets all zeros, 'read' requires --read_pqr", default=None, ) + config_group.add_argument( + "--radical_cofactor", + help="if the receptor has a radical cofactor and meeko faisl, try this option", + action="store_true", + ) box_group = parser.add_argument_group("Size and center of grid box") box_group.add_argument( @@ -590,6 +594,8 @@ def main(): mk_config["compute_charges"] = args.compute_charges + mk_config["radical_cofactor"] = args.radical_cofactor + # update config by inputs from arguments if args.charge_model is not None: diff --git a/meeko/polymer.py b/meeko/polymer.py index 9aaf99e6..8cd84591 100644 --- a/meeko/polymer.py +++ b/meeko/polymer.py @@ -1243,7 +1243,9 @@ def __init__( if unbound_unknown_res: for resname in set(unbound_unknown_res.values()): try: - cc = build_noncovalent_CC(resname, ignore_https_cert=ignore_https_cert) + cc = build_noncovalent_CC(resname, + ignore_https_cert=ignore_https_cert, + radical_cofactor=mk_prep.radical_cofactor) fetch_template_dict = json.loads(export_chem_templates_to_json([cc]))['residue_templates'][cc.resname] residue_templates.update({resname: ResidueTemplate( smiles = fetch_template_dict['smiles'], diff --git a/meeko/preparation.py b/meeko/preparation.py index fd87f2e1..443ec0d2 100644 --- a/meeko/preparation.py +++ b/meeko/preparation.py @@ -115,6 +115,7 @@ def __init__( crippen_as_solpar=False, override_ad4sol_par_including_q=False, override_ad4sol_par_including_q_qasp=0.0, + radical_cofactor=False, ): """ @@ -246,6 +247,8 @@ def __init__( self.add_index_map = add_index_map self.remove_smiles = remove_smiles + self.radical_cofactor=radical_cofactor + self._bond_typer = BondTyperLegacy() self._macrocycle_typer = FlexMacrocycle( self.min_ring_size, diff --git a/test/chemtempgen_test.py b/test/chemtempgen_test.py index 27618c2a..b31648a1 100644 --- a/test/chemtempgen_test.py +++ b/test/chemtempgen_test.py @@ -55,6 +55,16 @@ def test_build_noncovalent_CC(): assert template_equality_check(default_template_file, basename, "_fl-ccd", cc) + +def test_radical_cofactor_NO(): + basename = "NO" # free ligand from CCD + cc = build_noncovalent_CC(basename, radical_cofactor=True) + + assert cc is not None + assert isinstance(cc, ChemicalComponent) + + assert cc.smiles_exh == '[N]=O' + def test_add_variants(): basename = "AMP" cc_list = build_linked_CCs(basename)