Skills › Scientific Agent Skills › bioservices › compound_cross_reference.py
compound_cross_reference.py
424 lines · skills/bioservices/scripts/compound_cross_reference.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Compound Cross-Database Search
5This script searches for a compound by name and retrieves identifiers
6from multiple databases:
7- KEGG Compound
8- ChEBI
9- ChEMBL (via UniChem)
10- Basic compound properties
12Usage:
13 python compound_cross_reference.py COMPOUND_NAME [--output FILE]
15Examples:
16 python compound_cross_reference.py Geldanamycin
17 python compound_cross_reference.py "Adenosine triphosphate"
18 python compound_cross_reference.py Aspirin --output aspirin_info.txt
19"""
21import sys
22import argparse
23from urllib.parse import quote
24from bioservices import KEGG, UniChem, ChEBI, ChEMBL
27def search_kegg_compound(compound_name):
28 """Search KEGG for compound by name."""
29 print(f"\n{'='*70}")
30 print("STEP 1: KEGG Compound Search")
31 print(f"{'='*70}")
33 k = KEGG()
34 k.services.url = "https://rest.kegg.jp"
36 print(f"Searching KEGG for: {compound_name}")
38 try:
39 results = k.services.http_get("find/compound/" + quote(compound_name, safe=""), frmt="txt")
41 if not isinstance(results, str):
42 raise ValueError(f"KEGG search failed: {results!r}")
43 if not results.strip():
44 print(f"[FAIL] No results found in KEGG")
45 return k, None
47 # Parse results
48 lines = results.strip().split("\n")
49 print(f"[OK] Found {len(lines)} result(s):\n")
51 for i, line in enumerate(lines[:5], 1):
52 parts = line.split("\t")
53 kegg_id = parts[0]
54 description = parts[1] if len(parts) > 1 else "No description"
55 print(f" {i}. {kegg_id}: {description}")
57 exact_matches = [line for line in lines if len(line.split("\t", 1)) == 2
58 and compound_name.casefold() in {
59 name.strip().casefold() for name in line.split("\t", 1)[1].split(";")}]
60 if len(exact_matches) == 1:
61 lines = exact_matches
62 if len(lines) > 1:
63 print("[SKIP] Ambiguous compound name; refine the query before mapping")
64 return k, None
66 # Use the unique result
67 first_result = lines[0].split("\t")
68 kegg_id = first_result[0].replace("cpd:", "")
70 print(f"\nUsing: {kegg_id}")
72 return k, kegg_id
74 except Exception as e:
75 print(f"[FAIL] Error: {e}")
76 return k, None
79def get_kegg_info(kegg, kegg_id):
80 """Retrieve detailed KEGG compound information."""
81 print(f"\n{'='*70}")
82 print("STEP 2: KEGG Compound Details")
83 print(f"{'='*70}")
85 try:
86 print(f"Retrieving KEGG entry for {kegg_id}...")
88 entry = kegg.get(f"cpd:{kegg_id}")
90 if not entry:
91 print("[FAIL] Failed to retrieve entry")
92 return None
94 # Parse entry
95 compound_info = {
96 'kegg_id': kegg_id,
97 'name': None,
98 'formula': None,
99 'exact_mass': None,
100 'mol_weight': None,
101 'chebi_id': None,
102 'chebi_ids': [],
103 'pathways': []
104 }
106 current_section = None
108 for line in entry.split("\n"):
109 # KEGG field names start in column 1 and their continuation lines are
110 # indented. A new field therefore ends any multi-line section --
111 # without this reset, the indented DBLINKS lines that follow PATHWAY
112 # would be collected as pathways.
113 if line and not line.startswith(" "):
114 current_section = None
116 if line.startswith("NAME"):
117 compound_info['name'] = line.replace("NAME", "").strip().rstrip(";")
119 elif line.startswith("FORMULA"):
120 compound_info['formula'] = line.replace("FORMULA", "").strip()
122 elif line.startswith("EXACT_MASS"):
123 compound_info['exact_mass'] = line.replace("EXACT_MASS", "").strip()
125 elif line.startswith("MOL_WEIGHT"):
126 compound_info['mol_weight'] = line.replace("MOL_WEIGHT", "").strip()
128 elif "ChEBI:" in line:
129 parts = line.split("ChEBI:")
130 if len(parts) > 1:
131 for identifier in parts[1].strip().split():
132 if identifier not in compound_info['chebi_ids']:
133 compound_info['chebi_ids'].append(identifier)
135 elif line.startswith("PATHWAY"):
136 current_section = "pathway"
137 pathway = line.replace("PATHWAY", "").strip()
138 if pathway:
139 compound_info['pathways'].append(pathway)
141 elif current_section == "pathway" and line.startswith(" "):
142 pathway = line.strip()
143 if pathway:
144 compound_info['pathways'].append(pathway)
146 elif line.startswith(" ") and not line.startswith(" "):
147 current_section = None
149 # A KEGG entry can link distinct protonation or structural forms.
150 # Preserve every candidate; the legacy singular field is unique only.
151 if len(compound_info['chebi_ids']) == 1:
152 compound_info['chebi_id'] = compound_info['chebi_ids'][0]
153 elif compound_info['chebi_ids']:
154 print("[SKIP] Multiple ChEBI cross-references; review structures before mapping: "
155 + ", ".join(compound_info['chebi_ids']))
157 # Display information
158 print(f"\n[OK] KEGG Compound Information:")
159 print(f" ID: {compound_info['kegg_id']}")
160 print(f" Name: {compound_info['name']}")
161 print(f" Formula: {compound_info['formula']}")
162 print(f" Exact Mass: {compound_info['exact_mass']}")
163 print(f" Molecular Weight: {compound_info['mol_weight']}")
165 if compound_info['chebi_id']:
166 print(f" ChEBI ID: {compound_info['chebi_id']}")
168 if compound_info['pathways']:
169 print(f" Pathways: {len(compound_info['pathways'])} found")
171 return compound_info
173 except Exception as e:
174 print(f"[FAIL] Error: {e}")
175 return None
178def get_chembl_id(kegg_id, chebi_id=None):
179 """Resolve a KEGG entry's ChEBI cross-reference with UniChem 2.
181 KEGG is not a supported UniChem 2 source. Preserve the original KEGG ID
182 for provenance and use only the ChEBI ID recorded in that KEGG entry.
183 """
184 print(f"\n{'='*70}")
185 print("STEP 3: ChEMBL Mapping (via UniChem)")
186 print(f"{'='*70}")
187 if not chebi_id:
188 print(f"[SKIP] KEGG:{kegg_id} has no unique ChEBI cross-reference; mapping unavailable")
189 return None
190 chebi_id = str(chebi_id)
191 if not chebi_id.startswith("CHEBI:"):
192 chebi_id = f"CHEBI:{chebi_id}"
193 try:
194 response = UniChem().get_compounds(chebi_id, "chebi")
195 identifiers = sorted({
196 source["compoundId"]
197 for match in response.get("compounds", [])
198 for source in match.get("sources", [])
199 if source.get("shortName") == "chembl" and source.get("compoundId")
200 })
201 if len(identifiers) == 1:
202 print(f"[OK] {chebi_id} -> ChEMBL ID: {identifiers[0]}")
203 return identifiers[0]
204 if identifiers:
205 print(f"[SKIP] Ambiguous ChEMBL mappings: {', '.join(identifiers)}; review structures")
206 else:
207 print("[SKIP] No ChEMBL mapping returned; this does not prove absence")
208 except Exception as error:
209 print(f"[FAIL] UniChem lookup failed: {error}")
210 return None
213def get_chebi_info(chebi_id):
214 """Retrieve ChEBI compound information."""
215 print(f"\n{'='*70}")
216 print("STEP 4: ChEBI Details")
217 print(f"{'='*70}")
219 if not chebi_id:
220 print("[SKIP] No ChEBI ID available")
221 return None
223 try:
224 c = ChEBI()
226 print(f"Retrieving ChEBI entry for {chebi_id}...")
228 # Ensure proper format
229 if not chebi_id.startswith("CHEBI:"):
230 chebi_id = f"CHEBI:{chebi_id}"
232 entity = c.getCompleteEntity(chebi_id)
234 if entity:
235 print(f"\n[OK] ChEBI Information:")
236 print(f" ID: {entity.chebiId}")
237 print(f" Name: {entity.chebiAsciiName}")
239 if hasattr(entity, 'formula') and entity.formula:
240 print(f" Formula: {entity.formula}")
242 if hasattr(entity, 'mass') and entity.mass:
243 print(f" Mass: {entity.mass}")
245 if hasattr(entity, 'charge') and entity.charge:
246 print(f" Charge: {entity.charge}")
248 return {
249 'chebi_id': entity.chebiId,
250 'name': entity.chebiAsciiName,
251 'formula': entity.formula if hasattr(entity, 'formula') else None,
252 'mass': entity.mass if hasattr(entity, 'mass') else None
253 }
254 else:
255 print("[FAIL] Failed to retrieve ChEBI entry")
256 return None
258 except Exception as e:
259 print(f"[FAIL] Error: {e}")
260 return None
263def get_chembl_info(chembl_id):
264 """Retrieve ChEMBL compound information."""
265 print(f"\n{'='*70}")
266 print("STEP 5: ChEMBL Details")
267 print(f"{'='*70}")
269 if not chembl_id:
270 print("[SKIP] No ChEMBL ID available")
271 return None
273 try:
274 c = ChEMBL()
276 print(f"Retrieving ChEMBL entry for {chembl_id}...")
278 # `get_compound_by_chemblId` was a pre-1.6 name; current releases expose
279 # the same lookup as `get_molecule`.
280 compound = c.get_molecule(chembl_id)
282 if compound:
283 print(f"\n[OK] ChEMBL Information:")
284 print(f" ID: {chembl_id}")
286 if 'pref_name' in compound and compound['pref_name']:
287 print(f" Preferred Name: {compound['pref_name']}")
289 if 'molecule_properties' in compound:
290 props = compound['molecule_properties'] or {}
292 if 'full_mwt' in props:
293 print(f" Molecular Weight: {props['full_mwt']}")
295 if 'alogp' in props:
296 print(f" LogP: {props['alogp']}")
298 if 'hba' in props:
299 print(f" H-Bond Acceptors: {props['hba']}")
301 if 'hbd' in props:
302 print(f" H-Bond Donors: {props['hbd']}")
304 if 'molecule_structures' in compound:
305 structs = compound['molecule_structures'] or {}
307 if isinstance(structs.get('canonical_smiles'), str):
308 smiles = structs['canonical_smiles']
309 print(f" SMILES: {smiles[:60]}{'...' if len(smiles) > 60 else ''}")
311 return compound
312 else:
313 print("[FAIL] Failed to retrieve ChEMBL entry")
314 return None
316 except Exception as e:
317 print(f"[FAIL] Error: {e}")
318 return None
321def save_results(compound_name, kegg_info, chembl_id, output_file):
322 """Save results to file."""
323 print(f"\n{'='*70}")
324 print(f"Saving results to {output_file}")
325 print(f"{'='*70}")
327 with open(output_file, 'w') as f:
328 f.write("=" * 70 + "\n")
329 f.write(f"Compound Cross-Reference Report: {compound_name}\n")
330 f.write("=" * 70 + "\n\n")
332 # KEGG information
333 if kegg_info:
334 f.write("KEGG Compound\n")
335 f.write("-" * 70 + "\n")
336 f.write(f"ID: {kegg_info['kegg_id']}\n")
337 f.write(f"Name: {kegg_info['name']}\n")
338 f.write(f"Formula: {kegg_info['formula']}\n")
339 f.write(f"Exact Mass: {kegg_info['exact_mass']}\n")
340 f.write(f"Molecular Weight: {kegg_info['mol_weight']}\n")
341 f.write(f"Pathways: {len(kegg_info['pathways'])} found\n")
342 f.write("\n")
344 # Database IDs
345 f.write("Cross-Database Identifiers\n")
346 f.write("-" * 70 + "\n")
347 if kegg_info:
348 f.write(f"KEGG: {kegg_info['kegg_id']}\n")
349 if kegg_info['chebi_id']:
350 f.write(f"ChEBI: {kegg_info['chebi_id']}\n")
351 elif kegg_info.get('chebi_ids'):
352 f.write("ChEBI candidates (unresolved): "
353 + ", ".join(kegg_info['chebi_ids']) + "\n")
354 if chembl_id:
355 f.write(f"ChEMBL: {chembl_id}\n")
356 f.write("\n")
358 print(f"[OK] Results saved")
361def main():
362 """Main workflow."""
363 parser = argparse.ArgumentParser(
364 description="Search compound across multiple databases",
365 formatter_class=argparse.RawDescriptionHelpFormatter,
366 epilog="""
367Examples:
368 python compound_cross_reference.py Geldanamycin
369 python compound_cross_reference.py "Adenosine triphosphate"
370 python compound_cross_reference.py Aspirin --output aspirin_info.txt
371 """
372 )
373 parser.add_argument("compound", help="Compound name to search")
374 parser.add_argument("--output", default=None,
375 help="Output file for results (optional)")
377 args = parser.parse_args()
379 print("=" * 70)
380 print("BIOSERVICES: Compound Cross-Database Search")
381 print("=" * 70)
383 # Step 1: Search KEGG
384 kegg, kegg_id = search_kegg_compound(args.compound)
385 if not kegg_id:
386 print("\n[FAIL] Failed to find compound. Exiting.")
387 sys.exit(1)
389 # Step 2: Get KEGG details
390 kegg_info = get_kegg_info(kegg, kegg_id)
392 # Step 3: Map to ChEMBL
393 chembl_id = get_chembl_id(kegg_id, (kegg_info or {}).get("chebi_id"))
395 # Step 4: Get ChEBI details
396 chebi_info = None
397 if kegg_info and kegg_info['chebi_id']:
398 chebi_info = get_chebi_info(kegg_info['chebi_id'])
400 # Step 5: Get ChEMBL details
401 chembl_info = None
402 if chembl_id:
403 chembl_info = get_chembl_info(chembl_id)
405 # Summary
406 print(f"\n{'='*70}")
407 print("SUMMARY")
408 print(f"{'='*70}")
409 print(f" Compound: {args.compound}")
410 if kegg_info:
411 print(f" KEGG ID: {kegg_info['kegg_id']}")
412 if kegg_info['chebi_id']:
413 print(f" ChEBI ID: {kegg_info['chebi_id']}")
414 if chembl_id:
415 print(f" ChEMBL ID: {chembl_id}")
416 print(f"{'='*70}")
418 # Save to file if requested
419 if args.output:
420 save_results(args.compound, kegg_info, chembl_id, args.output)
423if __name__ == "__main__":
424 main()