#!/usr/bin/env python3 """ Generate ChEBI bulk submission SDF file from filtered glycan data. Combines glycan metadata with PubChem chemical structures. """ import pandas as pd import os import sys from datetime import datetime # PRODUCTION MODE today = datetime.now().strftime("%Y_%m_%d") data_dir = f"/data/projects/glygen/downloads/chebi_bulk_sub/{today}" input_csv = os.path.join(data_dir, "glycans_for_chebi.csv") pubchem_compounds_sdf = os.path.join(data_dir, "pubchem_glytoucan_compounds.sdf") pubchem_substances_sdf = os.path.join(data_dir, "pubchem_glytoucan_substances.sdf") output_sdf = os.path.join(data_dir, "glycans_chebi_bulk_submission.sdf") print(f"Input CSV: {input_csv}") print(f"PubChem compounds SDF: {pubchem_compounds_sdf}") print(f"PubChem substances SDF: {pubchem_substances_sdf}") print(f"Output SDF: {output_sdf}") print() # Check input files exist if not os.path.exists(input_csv): raise FileNotFoundError(f"ERROR: Input CSV not found: {input_csv}") if not os.path.exists(pubchem_compounds_sdf): raise FileNotFoundError(f"ERROR: PubChem compounds SDF not found: {pubchem_compounds_sdf}") if not os.path.exists(pubchem_substances_sdf): raise FileNotFoundError(f"ERROR: PubChem substances SDF not found: {pubchem_substances_sdf}") def parse_pubchem_sdf(sdf_file): """ Parse PubChem SDF file and create a dictionary mapping CID to structure block. Args: sdf_file: Path to the PubChem SDF file Returns: dict: {cid: structure_block_text} """ print(f"Parsing {os.path.basename(sdf_file)}...") structures = {} current_cid = None current_block = [] in_structure = False with open(sdf_file, 'r', encoding='utf-8', errors='ignore') as f: for line in f: # Start of new record (after $$$$) if line.startswith('$$$$'): if current_cid and current_block: # Save the structure block (everything before first >) structures[current_cid] = ''.join(current_block) current_cid = None current_block = [] in_structure = False continue # Look for PUBCHEM_COMPOUND_CID or PUBCHEM_SUBSTANCE_SID if line.startswith('> ') or line.startswith('> '): in_structure = False continue # If we see a CID/SID value line (right after the field name) if current_cid is None and line.strip() and not line.startswith('>'): # Check if previous lines indicate we're reading a CID/SID if len(current_block) == 0: # This might be the CID/SID - try to parse it try: potential_cid = line.strip() if potential_cid.isdigit(): current_cid = potential_cid current_block = [] in_structure = True except: pass # Collect structure block lines (before first >) if in_structure and not line.startswith('>'): current_block.append(line) elif line.startswith('>'): in_structure = False print(f" Found {len(structures)} structures") return structures def get_structure_block(pubchem_id, pubchem_type, compounds_dict, substances_dict): """ Get structure block for a given PubChem ID. Removes the first line (CID/SID number) but keeps the blank line to preserve structure. Args: pubchem_id: PubChem CID or SID pubchem_type: 'glycan_xref_pubchem_compound' or 'glycan_xref_pubchem_substance' compounds_dict: Dictionary of compound structures substances_dict: Dictionary of substance structures Returns: str: Structure block text with CID line removed, or empty block if not found """ pubchem_id_str = str(int(pubchem_id)) # Ensure it's a string without decimals structure_block = None if pubchem_type == 'glycan_xref_pubchem_compound': if pubchem_id_str in compounds_dict: structure_block = compounds_dict[pubchem_id_str] elif pubchem_type == 'glycan_xref_pubchem_substance': if pubchem_id_str in substances_dict: structure_block = substances_dict[pubchem_id_str] # If structure block found, remove the first line (CID/SID) but keep it blank if structure_block: lines = structure_block.split('\n') if lines: # Replace the first line (CID number) with a blank line lines[0] = '' structure_block = '\n'.join(lines) return structure_block # Return empty structure block if not found (starts with blank line for consistency) return "\n -OEChem-01312510302D\n\n 0 0 0 0 0 0 0 0 0999 V2000\nM END\n" def generate_sdf_entry(row, compounds_dict, substances_dict): """ Generate a single SDF entry for a glycan. Args: row: DataFrame row containing glycan data compounds_dict: Dictionary of compound structures substances_dict: Dictionary of substance structures Returns: str: Complete SDF entry """ # Get structure block structure_block = get_structure_block( row['pubchem_id'], row['pubchem_type'], compounds_dict, substances_dict ) # Build SDF entry sdf_lines = [] # Structure block sdf_lines.append(structure_block) # ID sdf_lines.append(f"> \n") sdf_lines.append(f"TEMP_{row['glytoucan_ac']}\n") sdf_lines.append("\n") # DEFINITION sdf_lines.append(f"> \n") sdf_lines.append(f"{row['definition']}\n") sdf_lines.append("\n") # NAME sdf_lines.append(f"> \n") sdf_lines.append(f"GlyTouCan {row['glytoucan_ac']}\n") sdf_lines.append("\n") # SYNONYM (WURCS) - only if available if pd.notna(row['sequence_wurcs']): sdf_lines.append(f"> \n") sdf_lines.append(f"{row['sequence_wurcs']}\n") sdf_lines.append("\n") # IUPAC_NAME - only if available if pd.notna(row['sequence_iupac_extended']): sdf_lines.append(f"> \n") sdf_lines.append(f"{row['sequence_iupac_extended']}\n") sdf_lines.append("\n") # RELATIONSHIP sdf_lines.append(f"> \n") sdf_lines.append(f"{row['relationship']}\n") sdf_lines.append("\n") # DATABASE_ACCESSION sdf_lines.append(f"> \n") accessions = f"GlyGen:{row['glytoucan_ac']}; GlyTouCan:{row['glytoucan_ac']}; PubChem:{int(row['pubchem_id'])}" sdf_lines.append(f"{accessions}\n") sdf_lines.append("\n") # REFERENCE (PubMed IDs) - only if available if pd.notna(row['pubmed_ids']): sdf_lines.append(f"> \n") # Convert pipe-separated to semicolon-space-separated pubmed_refs = row['pubmed_ids'].replace('|', '; ') sdf_lines.append(f"{pubmed_refs}\n") sdf_lines.append("\n") # COMMENT - currently empty as no DOIs available # Uncomment if needed in the future # sdf_lines.append(f"> \n") # sdf_lines.append(f"\n") # sdf_lines.append("\n") # End of record sdf_lines.append("$$$$\n") return ''.join(sdf_lines) # Main processing print("Reading glycan data...") df = pd.read_csv(input_csv) print(f"Loaded {len(df)} glycans\n") # Parse PubChem SDF files compounds_dict = parse_pubchem_sdf(pubchem_compounds_sdf) substances_dict = parse_pubchem_sdf(pubchem_substances_sdf) print() # Generate SDF file print(f"Generating SDF file: {output_sdf}") with open(output_sdf, 'w', encoding='utf-8') as f: for idx, row in df.iterrows(): if (idx + 1) % 100 == 0: print(f" Processed {idx + 1}/{len(df)} glycans...") sdf_entry = generate_sdf_entry(row, compounds_dict, substances_dict) f.write(sdf_entry) print(f"\nDone! Generated SDF with {len(df)} entries") print(f"Output file: {output_sdf}") # Show statistics structures_found = 0 structures_missing = 0 for idx, row in df.iterrows(): pubchem_id_str = str(int(row['pubchem_id'])) if row['pubchem_type'] == 'glycan_xref_pubchem_compound': if pubchem_id_str in compounds_dict: structures_found += 1 else: structures_missing += 1 else: if pubchem_id_str in substances_dict: structures_found += 1 else: structures_missing += 1 print(f"\nStructure statistics:") print(f" Structures found: {structures_found}") print(f" Structures missing (empty blocks): {structures_missing}")