Skills › Scientific Agent Skills › rdkit › substructure_filter.py
substructure_filter.py
380 lines · skills/rdkit/scripts/substructure_filter.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Substructure Filter
5Filter molecules based on substructure patterns using SMARTS.
6Supports inclusion and exclusion filters, and custom pattern libraries.
8Usage:
9 python substructure_filter.py molecules.smi --pattern "c1ccccc1" --output filtered.smi
10 python substructure_filter.py database.sdf --exclude "C(=O)Cl" --filter-type functional-groups
11"""
13import argparse
14import sys
15from pathlib import Path
17try:
18 from rdkit import Chem
19except ImportError:
20 print("Error: RDKit not installed. Install with: uv pip install rdkit")
21 sys.exit(1)
24from _common import SOURCE_INDEX, read_molecule_records, molecule_smiles
27# Neutral functional-group motifs; not universal chemistry or toxicity classifiers.
28PATTERN_LIBRARIES = {
29 'functional-groups': {
30 'alcohol': '[OX2H1][CX4]',
31 'aldehyde': '[CX3H1](=O)[#6]',
32 'ketone': '[#6][CX3](=O)[#6]',
33 'carboxylic_acid': 'C(=O)[OH]',
34 'ester': '[CX3](=O)[OX2H0][#6]',
35 'amide': 'C(=O)N',
36 'amine': '[NX3;$(N-[#6]);!$(N-C=O);!$(N-S(=O)=O);!$(N-P=O);!$(N=*)]',
37 'ether': '[#6][OX2H0;!$(O-C=O)][#6]',
38 'nitrile': 'C#N',
39 'nitro': '[N+](=O)[O-]',
40 'halide': '[CX4][F,Cl,Br,I]',
41 'thiol': '[#6][SX2H1]',
42 'sulfide': '[#6][SX2H0][#6]',
43 },
44 'rings': {
45 'benzene': 'c1ccccc1',
46 'pyridine': 'n1ccccc1',
47 'pyrrole': 'n1cccc1',
48 'furan': 'o1cccc1',
49 'thiophene': 's1cccc1',
50 'imidazole': 'n1cncc1',
51 'indole': 'c1ccc2[nH]ccc2c1',
52 'naphthalene': 'c1ccc2ccccc2c1',
53 },
54 'pains': {
55 'rhodanine': 'O=C1CSC(=S)N1',
56 'catechol': 'c1ccc(O)c(O)c1',
57 'quinone': 'O=C1C=CC(=O)C=C1',
58 'michael_acceptor': 'C=CC(=O)',
59 'alkyl_halide': '[C][I,Br]',
60 },
61 'privileged': {
62 'biphenyl': 'c1ccccc1-c2ccccc2',
63 'piperazine': 'N1CCNCC1',
64 'piperidine': 'N1CCCCC1',
65 'morpholine': 'N1CCOCC1',
66 }
67}
70def load_molecules(file_path, keep_props=True):
71 """Load molecules from file."""
72 # keep_props retained for callers; RDKit reads SDF properties by default.
73 return [mol for _, mol in read_molecule_records(file_path)]
76def create_pattern_query(pattern_string):
77 """Create SMARTS query from string or SMILES."""
78 # Try as SMARTS first
79 query = Chem.MolFromSmarts(pattern_string)
80 if query is not None and query.GetNumAtoms() > 0:
81 return query
83 # Try as SMILES
84 query = Chem.MolFromSmiles(pattern_string)
85 if query is not None and query.GetNumAtoms() > 0:
86 return query
88 print(f"Error: Invalid pattern: {pattern_string}")
89 return None
92def filter_molecules(molecules, include_patterns=None, exclude_patterns=None,
93 match_all_include=False, use_chirality=False):
94 """
95 Filter molecules based on substructure patterns.
97 Args:
98 molecules: List of RDKit Mol objects
99 include_patterns: List of (name, pattern) tuples to include
100 exclude_patterns: List of (name, pattern) tuples to exclude
101 match_all_include: If True, molecule must match ALL include patterns
103 Returns:
104 Tuple of (filtered_molecules, match_info)
105 """
106 for patterns in (include_patterns, exclude_patterns):
107 for name, pattern in patterns or []:
108 if pattern is None or pattern.GetNumAtoms() == 0:
109 raise ValueError(f"Invalid or empty query: {name}")
110 filtered = []
111 match_info = []
113 for idx, mol in enumerate(molecules):
114 if mol is None:
115 continue
117 source_index = mol.GetIntProp(SOURCE_INDEX) if mol.HasProp(SOURCE_INDEX) else idx + 1
119 # Check exclusion patterns first
120 excluded = False
121 exclude_matches = []
122 if exclude_patterns:
123 for name, pattern in exclude_patterns:
124 if mol.HasSubstructMatch(pattern, useChirality=use_chirality):
125 excluded = True
126 exclude_matches.append(name)
128 if excluded:
129 match_info.append({
130 'index': source_index,
131 'smiles': molecule_smiles(mol),
132 'status': 'excluded',
133 'matches': exclude_matches
134 })
135 continue
137 # Check inclusion patterns
138 if include_patterns:
139 include_matches = []
140 for name, pattern in include_patterns:
141 if mol.HasSubstructMatch(pattern, useChirality=use_chirality):
142 include_matches.append(name)
144 # Decide if molecule passes inclusion filter
145 if match_all_include:
146 passed = len(include_matches) == len(include_patterns)
147 else:
148 passed = len(include_matches) > 0
150 if passed:
151 filtered.append(mol)
152 match_info.append({
153 'index': source_index,
154 'smiles': molecule_smiles(mol),
155 'status': 'included',
156 'matches': include_matches
157 })
158 else:
159 match_info.append({
160 'index': source_index,
161 'smiles': molecule_smiles(mol),
162 'status': 'no_match',
163 'matches': []
164 })
165 else:
166 # No inclusion patterns, keep all non-excluded
167 filtered.append(mol)
168 match_info.append({
169 'index': source_index,
170 'smiles': molecule_smiles(mol),
171 'status': 'included',
172 'matches': []
173 })
175 return filtered, match_info
178def write_molecules(molecules, output_file):
179 """Write molecules to file."""
180 output_path = Path(output_file)
182 if output_path.suffix.lower() in ['.sdf']:
183 writer = Chem.SDWriter(str(output_path))
184 for mol in molecules:
185 writer.write(mol)
186 writer.close()
187 elif output_path.suffix.lower() in ['.smi', '.smiles', '.txt']:
188 with open(output_path, 'w') as f:
189 for mol in molecules:
190 smiles = molecule_smiles(mol)
191 name = mol.GetProp('_Name') if mol.HasProp('_Name') else ''
192 f.write(f"{smiles} {name}\n")
193 else:
194 raise ValueError(f'Unsupported output format: {output_path.suffix}')
196 print(f"Wrote {len(molecules)} molecules to {output_file}")
199def write_report(match_info, output_file):
200 """Write detailed match report."""
201 import csv
203 with open(output_file, 'w', newline='') as f:
204 fieldnames = ['Index', 'SMILES', 'Status', 'Matches']
205 writer = csv.DictWriter(f, fieldnames=fieldnames)
206 writer.writeheader()
208 for info in match_info:
209 writer.writerow({
210 'Index': info['index'],
211 'SMILES': info['smiles'],
212 'Status': info['status'],
213 'Matches': ', '.join(info['matches'])
214 })
217def print_summary(total, filtered, match_info):
218 """Print filtering summary."""
219 print("\n" + "="*60)
220 print("Filtering Summary")
221 print("="*60)
222 print(f"Total molecules: {total}")
223 print(f"Passed filter: {len(filtered)}")
224 print(f"Filtered out: {total - len(filtered)}")
225 print(f"Pass rate: {len(filtered)/total*100:.1f}%")
227 # Count by status
228 status_counts = {}
229 for info in match_info:
230 status = info['status']
231 status_counts[status] = status_counts.get(status, 0) + 1
233 print("\nBreakdown:")
234 for status, count in status_counts.items():
235 print(f" {status:15s}: {count}")
237 print("="*60)
240def main():
241 parser = argparse.ArgumentParser(
242 description='Filter molecules by substructure patterns',
243 formatter_class=argparse.RawDescriptionHelpFormatter,
244 epilog=f"""
245Pattern libraries:
246 --filter-type functional-groups Common functional groups
247 --filter-type rings Ring systems
248 --filter-type pains Five illustrative alerts, NOT the RDKit PAINS catalogue
249 --filter-type privileged Privileged structures
251Examples:
252 # Include molecules with benzene ring
253 python substructure_filter.py molecules.smi --pattern "c1ccccc1" -o filtered.smi
255 # Exclude reactive groups
256 python substructure_filter.py database.sdf --exclude "C(=O)Cl" -o clean.sdf
258 # Filter by functional groups
259 python substructure_filter.py molecules.smi --filter-type functional-groups -o fg.smi
261 # Exclude illustrative alert motifs (not a toxicity or PAINS verdict)
262 python substructure_filter.py compounds.smi --filter-type pains --exclude-mode -o clean.smi
264 # Multiple patterns
265 python substructure_filter.py mol.smi --pattern "c1ccccc1" --pattern "N" --match-all -o aromatic_amines.smi
266 """
267 )
269 parser.add_argument('input', nargs='?', help='Input file (SDF or SMILES)')
270 parser.add_argument('--pattern', '-p', action='append',
271 help='SMARTS/SMILES pattern to include (can specify multiple)')
272 parser.add_argument('--exclude', '-e', action='append',
273 help='SMARTS/SMILES pattern to exclude (can specify multiple)')
274 parser.add_argument('--filter-type', choices=PATTERN_LIBRARIES.keys(),
275 help='Use predefined pattern library')
276 parser.add_argument('--exclude-mode', action='store_true',
277 help='Use filter-type patterns for exclusion instead of inclusion')
278 parser.add_argument('--match-all', action='store_true',
279 help='Molecule must match ALL include patterns')
280 parser.add_argument('--chirality', action='store_true',
281 help='Require specified query stereochemistry')
282 parser.add_argument('--output', '-o', help='Output file')
283 parser.add_argument('--report', '-r', help='Write detailed report to CSV')
284 parser.add_argument('--list-patterns', action='store_true',
285 help='List available pattern libraries and exit')
287 args = parser.parse_args()
289 # List patterns mode
290 if args.list_patterns:
291 print("\nAvailable Pattern Libraries:")
292 print("="*60)
293 for lib_name, patterns in PATTERN_LIBRARIES.items():
294 print(f"\n{lib_name}:")
295 for name, pattern in patterns.items():
296 print(f" {name:25s}: {pattern}")
297 sys.exit(0)
299 if not args.input:
300 parser.error("input is required unless --list-patterns is used")
302 # Load molecules
303 print(f"Loading molecules from: {args.input}")
304 molecules = load_molecules(args.input)
305 if not molecules:
306 print("Error: No valid molecules loaded")
307 sys.exit(1)
309 print(f"Loaded {len(molecules)} molecules")
311 # Prepare patterns
312 include_patterns = []
313 exclude_patterns = []
315 # Add custom include patterns
316 if args.pattern:
317 for pattern_str in args.pattern:
318 query = create_pattern_query(pattern_str)
319 if query is None:
320 parser.error(f'Invalid include pattern: {pattern_str}')
321 include_patterns.append(('custom', query))
323 # Add custom exclude patterns
324 if args.exclude:
325 for pattern_str in args.exclude:
326 query = create_pattern_query(pattern_str)
327 if query is None:
328 parser.error(f'Invalid exclude pattern: {pattern_str}')
329 exclude_patterns.append(('custom', query))
331 # Add library patterns
332 if args.filter_type:
333 lib_patterns = PATTERN_LIBRARIES[args.filter_type]
334 for name, pattern_str in lib_patterns.items():
335 query = create_pattern_query(pattern_str)
336 if query:
337 if args.exclude_mode:
338 exclude_patterns.append((name, query))
339 else:
340 include_patterns.append((name, query))
342 if not include_patterns and not exclude_patterns:
343 print("Error: No patterns specified")
344 sys.exit(1)
346 # Print filter configuration
347 print(f"\nFilter configuration:")
348 if include_patterns:
349 print(f" Include patterns: {len(include_patterns)}")
350 if args.match_all:
351 print(" Mode: Match ALL")
352 else:
353 print(" Mode: Match ANY")
354 if exclude_patterns:
355 print(f" Exclude patterns: {len(exclude_patterns)}")
357 # Perform filtering
358 print("\nFiltering...")
359 filtered, match_info = filter_molecules(
360 molecules,
361 include_patterns=include_patterns if include_patterns else None,
362 exclude_patterns=exclude_patterns if exclude_patterns else None,
363 match_all_include=args.match_all,
364 use_chirality=args.chirality
365 )
367 # Print summary
368 print_summary(len(molecules), filtered, match_info)
370 # Write output
371 if args.output:
372 write_molecules(filtered, args.output)
374 if args.report:
375 write_report(match_info, args.report)
376 print(f"Detailed report written to: {args.report}")
379if __name__ == '__main__':
380 main()