#!/usr/bin/env python3 """ Post-process ChEBI SDF file to remove hydrogens. Uses the same method as ChEBI's libRDChEBI for consistency. """ import os import sys from chembl_structure_pipeline.standardizer import parse_molblock, update_mol_valences from rdkit import Chem # TESTING MODE - Read files from same directory as script # Change testing_mode to False to use production paths testing_mode = False if testing_mode: print("="*60) print("RUNNING IN TESTING MODE") print("Reading files from script directory") print("="*60) script_dir = os.path.dirname(os.path.abspath(__file__)) input_sdf = os.path.join(script_dir, "glycans_chebi_bulk_submission.sdf") output_sdf = os.path.join(script_dir, "glycans_chebi_bulk_submission_cleaned.sdf") else: # PRODUCTION MODE data_dir = f"/data/projects/glygen/downloads/chebi_bulk_sub/current/" input_sdf = os.path.join(data_dir, "glycans_chebi_bulk_submission.sdf") output_sdf = os.path.join(data_dir, "glycans_chebi_bulk_submission_cleaned.sdf") print(f"Input SDF: {input_sdf}") print(f"Output SDF: {output_sdf}") print() # Check input file exists if not os.path.exists(input_sdf): raise FileNotFoundError(f"ERROR: Input SDF not found: {input_sdf}") def transform_alias_to_r(molfile: str) -> str: """ Convert Carbon atoms with R-group aliases to proper R-group representations. Some molecules in old ChEBI have R groups defined as Carbons with aliases. This function converts them to proper R-group representations. Args: molfile (str): The molecule structure in molfile format Returns: str: Modified molecule structure in molfile format """ mol = parse_molblock(molfile) for at in mol.GetAtoms(): if "molFileAlias" in at.GetPropNames() and at.GetSymbol() == "C": alias = at.GetProp("molFileAlias") if alias.startswith("R"): at.SetAtomicNum(0) at.SetProp("dummyLabel", alias) at.SetProp("molFileAlias", "") return Chem.MolToMolBlock(mol) def remove_hs(molfile: str) -> str: """ Remove hydrogen atoms from a molecule structure. Bespoke remove Hs function for MetaboLights team. Preserves stereochemistry-relevant hydrogen atoms. Args: molfile (str): The molecule structure in molfile format Returns: str: Modified molecule structure in molfile format with hydrogens removed """ mol = parse_molblock(molfile) Chem.FastFindRings(mol) mol = update_mol_valences(mol) indices = [] for atom in mol.GetAtoms(): if atom.GetAtomicNum() == 1 and not atom.GetIsotope(): bnd = atom.GetBonds()[0] if ( bnd.GetBondDir() not in (Chem.BondDir.BEGINWEDGE, Chem.BondDir.BEGINDASH) ) and not ( bnd.HasProp("_MolFileBondStereo") and bnd.GetUnsignedProp("_MolFileBondStereo") in (1, 6) ): indices.append(atom.GetIdx()) mol = Chem.RWMol(mol) for index in sorted(indices, reverse=True): mol.RemoveAtom(index) props = molfile.split("M END")[1].strip() props = props if len(props) > 1 else None out_molfile = Chem.MolToMolBlock(mol) if props: out_molfile += props return out_molfile def process_sdf_entry(entry_lines): """ Process a single SDF entry to remove hydrogens. Args: entry_lines: List of lines for one SDF entry Returns: str: Processed SDF entry """ # Find where the mol block ends (M END line) mol_end_idx = None for i, line in enumerate(entry_lines): if line.strip() == 'M END': mol_end_idx = i break if mol_end_idx is None: # No mol block found (empty structure), return as-is return ''.join(entry_lines) # Extract mol block and metadata mol_block = ''.join(entry_lines[:mol_end_idx + 1]) metadata = ''.join(entry_lines[mol_end_idx + 1:]) # Check if this is an empty structure block (0 atoms) lines = mol_block.split('\n') if len(lines) >= 4: counts_line = lines[3].strip() # Parse atom count (first number in the counts line) parts = counts_line.split() if parts and parts[0] == '0': # Empty structure, return as-is return ''.join(entry_lines) try: # Transform R-group aliases mol_block = transform_alias_to_r(mol_block) # Remove hydrogens mol_block = remove_hs(mol_block) # Combine processed mol block with metadata return mol_block + metadata except Exception as e: print(f" WARNING: Failed to process structure: {e}") print(f" Keeping original structure") return ''.join(entry_lines) def process_sdf_file(input_path, output_path): """ Process entire SDF file to remove hydrogens from all entries. Args: input_path: Path to input SDF file output_path: Path to output SDF file """ print("Processing SDF file...") current_entry = [] entries_processed = 0 entries_with_structures = 0 entries_modified = 0 with open(input_path, 'r', encoding='utf-8', errors='ignore') as infile: with open(output_path, 'w', encoding='utf-8') as outfile: for line in infile: current_entry.append(line) # End of entry marker if line.strip() == '$$$$': entries_processed += 1 if entries_processed % 100 == 0: print(f" Processed {entries_processed} entries...") # Process the entry original_entry = ''.join(current_entry) processed_entry = process_sdf_entry(current_entry) # Check if structure was modified if original_entry != processed_entry: entries_modified += 1 # Check if it has a real structure (not empty) if ' 0 0 0 0 0 0 0 0 0999 V2000' not in processed_entry: entries_with_structures += 1 # Write processed entry outfile.write(processed_entry) # Reset for next entry current_entry = [] print(f"\nProcessing complete!") print(f" Total entries: {entries_processed}") print(f" Entries with structures: {entries_with_structures}") print(f" Entries modified (H removed): {entries_modified}") print(f" Empty structures: {entries_processed - entries_with_structures}") # Main processing try: process_sdf_file(input_sdf, output_sdf) print(f"\nOutput file: {output_sdf}") print("\nDone!") except Exception as e: print(f"\nERROR: {e}", file=sys.stderr) sys.exit(1)