{
  "skill": "alterlab-brenda",
  "evals": [
    {
      "id": "km-values-by-ec",
      "prompt": "I'm building a kinetic model of alcohol dehydrogenase (EC 1.1.1.1). Pull the Km values for ethanol across organisms from BRENDA so I can compare turnover.",
      "expected_output": "Invokes alterlab-brenda: authenticates against the BRENDA SOAP API and calls get_km_values(\"1.1.1.1\", substrate=\"ethanol\"), then parses the returned organism*/substrate*/kmValue* entries (parse_km_entry) into per-organism Km values with pH/temperature commentary.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "behavior", "value": "Uses the EC number 1.1.1.1 and retrieves Km kinetic values via the BRENDA SOAP API rather than a generic web search." }
      ]
    },
    {
      "id": "cross-organism-comparison",
      "prompt": "Compare the kinetic properties and optimal pH/temperature of glutamine synthetase across E. coli, S. cerevisiae, and Thermus thermophilus for an enzyme-engineering project.",
      "expected_output": "Invokes alterlab-brenda: uses compare_across_organisms with the relevant EC number plus get_environmental_parameters to report average Km, optimal pH, and temperature range per organism, drawing all data from BRENDA.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "behavior", "value": "Compares the same enzyme across multiple named organisms using BRENDA organism-specific kinetic and environmental data." }
      ]
    },
    {
      "id": "inhibitor-lookup",
      "prompt": "What inhibitors and Ki values are reported for lactate dehydrogenase in BRENDA? I need the regulation data for a drug-target writeup.",
      "expected_output": "Invokes alterlab-brenda: calls get_inhibitors on the EC number for lactate dehydrogenase and returns inhibitor names, inhibition type, Ki and IC50 values from BRENDA's enzyme regulation data.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "Ki" }
      ]
    },
    {
      "id": "thermophilic-homologs",
      "prompt": "I need a heat-stable version of an enzyme acting on ethanol for an industrial process. Find thermophilic homologs of EC 1.1.1.1 with optimal temperatures above 60 C.",
      "expected_output": "Invokes alterlab-brenda: uses find_thermophilic_homologs(\"1.1.1.1\", min_temp=60) and reports candidate organisms with their optimal temperatures and Km values for industrial enzyme selection.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "behavior", "value": "Filters BRENDA enzyme variants by optimal temperature to surface thermostable candidates." }
      ]
    },
    {
      "id": "enzyme-by-substrate-discovery",
      "prompt": "Which enzymes can act on 2-phenylethanol as a substrate? I want EC numbers and reaction equations to plan a biosynthetic route.",
      "expected_output": "Invokes alterlab-brenda: discovers candidate enzymes for 2-phenylethanol via BRENDA reaction data (search_by_pattern / get_reactions), whose entries carry ecNumber, and returns EC numbers and reaction equations. (Km-based search_enzymes_by_substrate alone cannot supply EC numbers, since getKmValue responses omit ecNumber.)",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "EC" }
      ]
    },
    {
      "id": "near-miss-kegg",
      "prompt": "Map out the full glycolysis pathway with its reaction steps and how the metabolites connect, so I can see the whole metabolic map.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-kegg. The user wants a curated metabolic pathway map and reaction network, not enzyme kinetic parameters (Km/kcat/Ki) for a specific EC number, which is BRENDA's domain.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-kegg" }
      ]
    },
    {
      "id": "near-miss-uniprot",
      "prompt": "Give me the amino acid sequence and domain architecture of human alcohol dehydrogenase ADH1B.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-uniprot. The user wants protein sequence and domain annotation, not enzyme kinetic measurements, which is what BRENDA provides.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-uniprot" }
      ]
    }
  ]
}

# BRENDA Database API Reference

## Overview

This document provides detailed reference information for the BRENDA (BRaunschweig ENzyme DAtabase) SOAP API and the Python client implementation. BRENDA is the world's most comprehensive enzyme information system, containing over 45,000 enzymes with millions of kinetic data points.

## SOAP API Endpoints

### Base WSDL URL
```
https://www.brenda-enzymes.org/soap/brenda_zeep.wsdl
```

### Authentication

All BRENDA API calls require authentication using email and password:

**Parameters:**
- `email`: Your registered BRENDA email address
- `password`: Your BRENDA account password

**Authentication Process:**
1. Password is hashed using SHA-256 before transmission
2. Email and hashed password are the first two arguments of every call
3. Legacy support for `BRENDA_EMIAL` environment variable (note the typo)

**Calling convention (zeep).** Pass each value as a separate positional argument in the
WSDL's part order (listed per action below), exactly as BRENDA's own example does:

```python
from zeep import Client, Settings
import hashlib

client = Client("https://www.brenda-enzymes.org/soap/brenda_zeep.wsdl",
                settings=Settings(strict=False))
pw = hashlib.sha256("myPassword".encode("utf-8")).hexdigest()
result = client.service.getKmValue("j.doe@example.edu", pw, "ecNumber*1.1.1.1",
                                   "organism*Homo sapiens", "kmValue*", "kmValueMaximum*",
                                   "substrate*", "commentary*", "ligandStructureId*",
                                   "literature*")
```

A single comma-joined string (the older SOAPpy style) raises
`ValidationError: Missing element password` in zeep. Inspect any operation's order with
`client.service._binding._operations["getKmValue"].input.signature()`.

## Available SOAP Actions

### getKmValue

Retrieves Michaelis constant (Km) values for enzymes.

**Parameters:**
1. `email`: BRENDA account email
2. `passwordHash`: SHA-256 hashed password
3. `ecNumber*: EC number of the enzyme (wildcards allowed)
4. `organism*: Organism name (wildcards allowed, default: "*")
5. `kmValue*: Km value field (default: "*")
6. `kmValueMaximum*: Maximum Km value field (default: "*")
7. `substrate*: Substrate name (wildcards allowed, default: "*")
8. `commentary*: Commentary field (default: "*")
9. `ligandStructureId*: Ligand structure ID field (default: "*")
10. `literature*: Literature reference field (default: "*")

**Wildcards:**
- `*`: Matches any sequence
- Can be used with partial EC numbers (e.g., "1.1.*")

**Response Format:** zeep returns a list of `kmValueObject` records (fields below;
`literature` is a list of BRENDA reference IDs). `brenda_client.split_entries()` renders
each as the legacy delimited string used by the parsers:
```
organism*Escherichia coli#substrate*glucose#kmValue*0.12#kmValueMaximum*#commentary*pH 7.4, 25°C#ligandStructureId*#literature*
```

**Example Response Fields:**
- `organism`: Source organism
- `substrate`: Substrate name
- `kmValue`: Michaelis constant value (typically in mM)
- `kmValueMaximum`: Maximum Km value (if available)
- `commentary`: Experimental conditions (pH, temperature, etc.)
- `ligandStructureId`: BRENDA ligand structure identifier
- `literature`: Reference to primary literature

### getReaction

Retrieves reaction equations and stoichiometry for enzymes.

**Parameters:**
1. `email`: BRENDA account email
2. `passwordHash`: SHA-256 hashed password
3. `ecNumber*: EC number of the enzyme (wildcards allowed)
4. `organism*: Organism name (wildcards allowed, default: "*")
5. `reaction*: Reaction equation (wildcards allowed, default: "*")
6. `commentary*: Commentary field (default: "*")
7. `literature*: Literature reference field (default: "*")

**Response Format:**
```
ecNumber*1.1.1.1#organism*Saccharomyces cerevisiae#reaction*ethanol + NAD+ <=> acetaldehyde + NADH + H+#commentary*#literature*
```

**Example Response Fields:**
- `ecNumber`: Enzyme Commission number
- `organism`: Source organism
- `reaction`: Balanced chemical equation (using <=> for equilibrium, -> for direction)
- `commentary`: Additional information
- `literature`: Reference citation

## Data Field Specifications

### EC Number Format

EC numbers follow the standard hierarchical format: `A.B.C.D`

- **A**: Main class (1-6)
  - 1: Oxidoreductases
  - 2: Transferases
  - 3: Hydrolases
  - 4: Lyases
  - 5: Isomerases
  - 6: Ligases
- **B**: Subclass
- **C**: Sub-subclass
- **D**: Serial number

**Examples:**
- `1.1.1.1`: Alcohol dehydrogenase
- `1.1.1.2`: Alcohol dehydrogenase (NADP+)
- `3.2.1.23`: Beta-galactosidase
- `2.7.1.1`: Hexokinase

### Organism Names

Organism names should use proper binomial nomenclature:

**Correct Format:**
- `Escherichia coli`
- `Saccharomyces cerevisiae`
- `Homo sapiens`

**Wildcards:**
- `Escherichia*`: Matches all E. coli strains
- `*coli`: Matches all coli species
- `*`: Matches all organisms

### Substrate Names

Substrate names follow IUPAC or common biochemical conventions:

**Common Formats:**
- Chemical names: `glucose`, `ethanol`, `pyruvate`
- IUPAC names: `β-D-glucose`, `ethanol`, `2-oxopropanoic acid`
- Abbreviations: `ATP`, `NAD+`, `CoA`

**Special Cases:**
- Cofactors: `NAD+`, `NADH`, `NADP+`, `NADPH`
- Metal ions: `Mg2+`, `Zn2+`, `Fe2+`
- Inorganic compounds: `H2O`, `CO2`, `O2`

### Commentary Field Format

Commentary fields contain experimental conditions and other metadata:

**Common Information:**
- **pH**: `pH 7.4`, `pH 6.5-8.0`
- **Temperature**: `25°C`, `37°C`, `50-60°C`
- **Buffer systems**: `phosphate buffer`, `Tris-HCl`
- **Purity**: `purified enzyme`, `crude extract`
- **Assay conditions**: `spectrophotometric`, `radioactive`
- **Inhibition**: `inhibited by heavy metals`, `activated by Mg2+`

**Examples:**
- `pH 7.4, 25°C, phosphate buffer`
- `pH 6.5-8.0 optimum, thermostable enzyme`
- `purified enzyme, specific activity 125 U/mg`
- `inhibited by iodoacetate, activated by Mn2+`

### Reaction Equation Format

Reactions use standard biochemical notation:

**Symbols:**
- `+`: Separate reactants/products
- `<=>`: Reversible reactions
- `->`: Irreversible (directional) reactions
- `=`: Alternative notation for reactions

**Common Patterns:**
- **Oxidation/reduction**: `alcohol + NAD+ <=> aldehyde + NADH + H+`
- **Phosphorylation**: `glucose + ATP <=> glucose-6-phosphate + ADP`
- **Hydrolysis**: `ester + H2O <=> acid + alcohol`
- **Carboxylation**: `acetyl-CoA + CO2 + H2O <=> malonyl-CoA`

**Cofactor Requirements:**
- **Oxidoreductases**: NAD+, NADH, NADP+, NADPH, FAD, FADH2
- **Transferases**: ATP, ADP, GTP, GDP
- **Ligases**: ATP, CoA

## Rate Limiting and Usage

### API Rate Limits

- BRENDA's SOAP page asks users not to send **more than one request per second**
  (faster clients may be treated as bots); `brenda_client.call_brenda()` enforces this
- Registration is required; data are licensed **CC BY 4.0**

### Best Practices

1. **Implement delays**: Add 0.5-1 second between requests
2. **Cache results**: Store frequently accessed data locally
3. **Use specific searches**: Narrow by organism and substrate when possible
4. **Batch operations**: Group related queries
5. **Handle errors gracefully**: Check for HTTP and SOAP errors
6. **Use wildcards judiciously**: Broad searches return large datasets

### Error Handling

**Common SOAP Errors:**
- `Authentication failed`: Check email/password
- `No data found`: Verify EC number, organism, substrate spelling
- `Rate limit exceeded`: Reduce request frequency
- `Invalid parameters`: Check parameter format and order

**Network Errors:**
- Connection timeouts
- SSL/TLS errors
- Service unavailable

## Python Client Reference

### brenda_client Module

#### Core Functions

**`load_env_from_file(path=".env")`**
- **Purpose**: Load environment variables from .env file
- **Parameters**: `path` - Path to .env file (default: ".env")
- **Returns**: None (populates os.environ)

**`_get_credentials() -> tuple[str, str]`**
- **Purpose**: Retrieve BRENDA credentials from environment
- **Returns**: Tuple of (email, password)
- **Raises**: RuntimeError if credentials missing

**`_get_client() -> Client`**
- **Purpose**: Initialize or retrieve SOAP client
- **Returns**: Zeep Client instance
- **Features**: Singleton pattern, custom transport settings

**`_hash_password(password: str) -> str`**
- **Purpose**: Generate SHA-256 hash of password
- **Parameters**: `password` - Plain text password
- **Returns**: Hexadecimal SHA-256 hash

**`call_brenda(action: str, parameters: List[str])`**
- **Purpose**: Execute BRENDA SOAP action (throttled to one call per second)
- **Parameters**:
  - `action` - SOAP action name (e.g., "getKmValue")
  - `parameters` - `field*value` tokens in the operation's WSDL order
- **Returns**: Raw zeep response (list of typed result objects); pass it through
  `split_entries()` for `field*value#…` strings

#### Convenience Functions

**`get_km_values(ec_number: str, organism: str = "*", substrate: str = "*") -> List[str]`**
- **Purpose**: Retrieve Km values for specified enzyme
- **Parameters**:
  - `ec_number`: Enzyme Commission number
  - `organism`: Organism name (wildcard allowed, default: "*")
  - `substrate`: Substrate name (wildcard allowed, default: "*")
- **Returns**: List of parsed data strings

**`get_reactions(ec_number: str, organism: str = "*", reaction: str = "*") -> List[str]`**
- **Purpose**: Retrieve reaction data for specified enzyme
- **Parameters**:
  - `ec_number`: Enzyme Commission number
  - `organism`: Organism name (wildcard allowed, default: "*")
  - `reaction`: Reaction pattern (wildcard allowed, default: "*")
- **Returns**: List of reaction data strings

#### Utility Functions

**`split_entries(return_text: str) -> List[str]`**
- **Purpose**: Normalize BRENDA responses to list format
- **Parameters**: `return_text` - Raw response from BRENDA
- **Returns**: List of individual data entries
- **Features**: Handles both string and complex object responses

## Data Structures and Parsing

### Km Entry Structure

**Parsed Km Entry Dictionary:**
```python
{
    'ecNumber': '1.1.1.1',
    'organism': 'Escherichia coli',
    'substrate': 'ethanol',
    'kmValue': '0.12',
    'km_value_numeric': 0.12,  # Extracted numeric value
    'kmValueMaximum': '',
    'commentary': 'pH 7.4, 25°C',
    'ph': 7.4,               # Extracted from commentary
    'temperature': 25.0,      # Extracted from commentary
    'ligandStructureId': '',
    'literature': ''
}
```

### Reaction Entry Structure

**Parsed Reaction Entry Dictionary:**
```python
{
    'ecNumber': '1.1.1.1',
    'organism': 'Saccharomyces cerevisiae',
    'reaction': 'ethanol + NAD+ <=> acetaldehyde + NADH + H+',
    'reactants': ['ethanol', 'NAD+'],
    'products': ['acetaldehyde', 'NADH', 'H+'],
    'commentary': '',
    'literature': ''
}
```

## Query Patterns and Examples

### Basic Queries

**Get all Km values for an enzyme:**
```python
from scripts.brenda_client import get_km_values

# Get all alcohol dehydrogenase Km values
km_data = get_km_values("1.1.1.1")
```

**Get Km values for specific organism:**
```python
# Get human alcohol dehydrogenase Km values
human_km = get_km_values("1.1.1.1", organism="Homo sapiens")
```

**Get Km values for specific substrate:**
```python
# Get Km for ethanol oxidation
ethanol_km = get_km_values("1.1.1.1", substrate="ethanol")
```

### Wildcard Searches

**Search for enzyme families:**
```python
# All alcohol dehydrogenases
alcohol_dehydrogenases = get_km_values("1.1.1.*")

# All hexokinases
hexokinases = get_km_values("2.7.1.*")
```

**Search for organism groups:**
```python
# All E. coli strains
e_coli_enzymes = get_km_values("*", organism="Escherichia coli")

# All Bacillus species
bacillus_enzymes = get_km_values("*", organism="Bacillus*")
```

### Combined Searches

**Specific enzyme-substrate combination:**
```python
# Get Km values for glucose oxidation in yeast
glucose_km = get_km_values("1.1.1.1",
                          organism="Saccharomyces cerevisiae",
                          substrate="glucose")
```

### Reaction Queries

**Get all reactions for an enzyme:**
```python
from scripts.brenda_client import get_reactions

reactions = get_reactions("1.1.1.1")
```

**Search for reactions with specific substrates:**
```python
# Find reactions involving glucose
glucose_reactions = get_reactions("*", reaction="*glucose*")
```

## Data Analysis Patterns

### Kinetic Parameter Analysis

**Extract numeric Km values:**
```python
from scripts.brenda_queries import parse_km_entry

km_data = get_km_values("1.1.1.1", substrate="ethanol")
numeric_kms = []

for entry in km_data:
    parsed = parse_km_entry(entry)
    if 'km_value_numeric' in parsed:
        numeric_kms.append(parsed['km_value_numeric'])

if numeric_kms:
    print(f"Average Km: {sum(numeric_kms)/len(numeric_kms):.3f}")
    print(f"Range: {min(numeric_kms):.3f} - {max(numeric_kms):.3f}")
```

### Organism Comparison

**Compare enzyme properties across organisms:**
```python
from scripts.brenda_queries import compare_across_organisms

organisms = ["Escherichia coli", "Saccharomyces cerevisiae", "Homo sapiens"]
comparison = compare_across_organisms("1.1.1.1", organisms)

for org_data in comparison:
    if org_data.get('data_points', 0) > 0:
        print(f"{org_data['organism']}: {org_data['average_km']:.3f}")
```

### Substrate Specificity

**Analyze substrate preferences:**
```python
from scripts.brenda_queries import get_substrate_specificity

specificity = get_substrate_specificity("1.1.1.1")

for substrate_data in specificity[:5]:  # Top 5
    print(f"{substrate_data['name']}: Km = {substrate_data['km']:.3f}")
```

## Integration Examples

### Metabolic Pathway Construction

**Build enzymatic pathway:**
```python
from scripts.enzyme_pathway_builder import find_pathway_for_product

# Find pathway for lactate production
pathway = find_pathway_for_product("lactate", max_steps=3)

for step in pathway['steps']:
    print(f"Step {step['step_number']}: {step['substrate']} -> {step['product']}")
    print(f"Enzymes available: {len(step['enzymes'])}")
```

### Enzyme Engineering Support

**Find thermostable variants:**
```python
from scripts.brenda_queries import find_thermophilic_homologs

thermophilic = find_thermophilic_homologs("1.1.1.1", min_temp=50)

for enzyme in thermophilic:
    print(f"{enzyme['organism']}: {enzyme['optimal_temperature']}°C")
```

### Kinetic Modeling

**Extract parameters for modeling:**
```python
from scripts.brenda_queries import get_modeling_parameters

model_data = get_modeling_parameters("1.1.1.1", substrate="ethanol")

print(f"Km: {model_data['km']}")
print(f"Vmax: {model_data['vmax']}")
print(f"Optimal conditions: pH {model_data['ph']}, {model_data['temperature']}°C")
```

## Troubleshooting

### Common Issues

**Authentication Errors:**
- Check BRENDA_EMAIL and BRENDA_PASSWORD environment variables
- Verify account is active and has API access
- Note legacy BRENDA_EMIAL support (typo in variable name)

**No Data Returned:**
- Verify EC number format (e.g., "1.1.1.1", not "1.1.1")
- Check spelling of organism and substrate names
- Try wildcards for broader searches
- Some enzymes may have limited data in BRENDA

**Rate Limiting:**
- Implement delays between requests
- Cache results locally
- Use more specific queries to reduce data volume
- Consider batch operations

**Data Format Issues:**
- Use provided parsing functions
- Handle missing fields gracefully
- BRENDA data format can be inconsistent
- Validate parsed data before use

### Performance Optimization

**Query Efficiency:**
- Use specific EC numbers when known
- Limit by organism or substrate to reduce result size
- Cache frequently accessed data
- Batch similar requests

**Memory Management:**
- Process large datasets in chunks
- Use generators for large result sets
- Clear parsed data when no longer needed

**Network Optimization:**
- Implement retry logic for network errors
- Use appropriate timeouts
- Monitor request frequency

## Additional Resources

### Official Documentation

- **BRENDA Website**: https://www.brenda-enzymes.org/
- **SOAP API Documentation**: https://www.brenda-enzymes.org/soap.php
- **Enzyme Nomenclature**: https://www.iubmb.org/enzyme/
- **EC Number Database**: https://www.qmul.ac.uk/sbcs/iubmb/enzyme/

### Related Libraries

- **Zeep (SOAP Client)**: https://python-zeep.readthedocs.io/
- **PubChemPy**: https://pubchempy.readthedocs.io/
- **BioPython**: https://biopython.org/
- **RDKit**: https://www.rdkit.org/

### Data Formats

- **Enzyme Commission Numbers**: IUBMB enzyme classification
- **IUPAC Nomenclature**: Chemical naming conventions
- **Biochemical Reactions**: Standard equation notation
- **Kinetic Parameters**: Michaelis-Menten kinetics

### Community Resources

- **BRENDA Help Desk**: Support via official website
- **Bioinformatics Forums**: Stack Overflow, Biostars
- **GitHub Issues**: Project-specific bug reports
- **Research Papers**: Primary literature for enzyme data

---

*This API reference covers the core functionality of the BRENDA SOAP API and Python client. For complete details on available data fields and query patterns, consult the official BRENDA documentation.*

# BRENDA — Capabilities & Code Examples

Copy-ready snippets for each BRENDA capability area. For the underlying SOAP API methods
and parameter lists, see `api_reference.md`. For end-to-end recipes, see `workflows.md`.
For the helper-script function inventory, see `helper_scripts.md`.

## 1. Kinetic Parameter Retrieval

**Get Km Values by EC Number**:
```python
from scripts.brenda_client import get_km_values

# Get Km values for all organisms
km_data = get_km_values("1.1.1.1")  # Alcohol dehydrogenase

# Get Km values for specific organism
km_data = get_km_values("1.1.1.1", organism="Saccharomyces cerevisiae")

# Get Km values for specific substrate
km_data = get_km_values("1.1.1.1", substrate="ethanol")
```

**Parse Km Results**:
```python
for entry in km_data:
    print(f"Km: {entry}")
    # Example output: "organism*Homo sapiens#substrate*ethanol#kmValue*1.2#commentary*"
```

**Extract Specific Information**:
```python
from scripts.brenda_queries import parse_km_entry, extract_organism_data

for entry in km_data:
    parsed = parse_km_entry(entry)
    organism = extract_organism_data(entry)
    # parse_km_entry exposes raw BRENDA field 'kmValue' (string) and the
    # derived 'km_value_numeric' (float); there is no 'km_value' key.
    print(f"Organism: {parsed.get('organism')}")
    print(f"Substrate: {parsed.get('substrate')}")
    print(f"Km value: {parsed.get('kmValue')}  (numeric: {parsed.get('km_value_numeric')})")
    print(f"pH: {parsed.get('ph', 'N/A')}")
    print(f"Temperature: {parsed.get('temperature', 'N/A')}")
```

## 2. Reaction Information

**Get Reactions by EC Number**:
```python
from scripts.brenda_client import get_reactions

# Get all reactions for EC number
reactions = get_reactions("1.1.1.1")

# Filter by organism
reactions = get_reactions("1.1.1.1", organism="Escherichia coli")

# Search specific reaction
reactions = get_reactions("1.1.1.1", reaction="ethanol + NAD+")
```

**Process Reaction Data**:
```python
from scripts.brenda_queries import parse_reaction_entry, extract_substrate_products

for reaction in reactions:
    parsed = parse_reaction_entry(reaction)
    substrates, products = extract_substrate_products(reaction)

    print(f"Reaction: {parsed['reaction']}")
    print(f"Organism: {parsed['organism']}")
    print(f"Substrates: {substrates}")
    print(f"Products: {products}")
```

## 3. Enzyme Discovery

> Gotcha: `getKmValue` responses do **not** carry an `ecNumber` field, so
> `search_enzymes_by_substrate` (built on Km data) returns `ec_number=''` and has
> no `enzyme_name`/`reaction` keys — it yields organism + substrate + Km. To resolve
> EC numbers and reactions for a substrate, query reaction data
> (`search_by_pattern` / `get_reactions`), whose entries include `ecNumber`.

**Find Enzymes by Substrate** (organism/substrate/Km, no EC number):
```python
from scripts.brenda_queries import search_enzymes_by_substrate

# Find enzymes that act on glucose
enzymes = search_enzymes_by_substrate("glucose", limit=20)

for enzyme in enzymes:
    print(f"Organism: {enzyme['organism']}")
    print(f"Substrate: {enzyme['substrate']}")
    print(f"Km: {enzyme['km_value']}")
```

**Find Enzymes by Product**:
```python
from scripts.brenda_queries import search_enzymes_by_product

# Find enzymes that produce lactate
enzymes = search_enzymes_by_product("lactate", limit=10)
```

**Search by Reaction Pattern**:
```python
from scripts.brenda_queries import search_by_pattern

# Find oxidation reactions
enzymes = search_by_pattern("oxidation", limit=15)
```

## 4. Organism-Specific Enzyme Data

**Get Enzyme Data for Multiple Organisms**:
```python
from scripts.brenda_queries import compare_across_organisms

organisms = ["Escherichia coli", "Saccharomyces cerevisiae", "Homo sapiens"]
comparison = compare_across_organisms("1.1.1.1", organisms)

for org_data in comparison:
    print(f"Organism: {org_data['organism']}")
    print(f"Avg Km: {org_data['average_km']}")
    print(f"Optimal pH: {org_data['optimal_ph']}")
    print(f"Temperature range: {org_data['temperature_range']}")
```

**Find Organisms with Specific Enzyme**:
```python
from scripts.brenda_queries import get_organisms_for_enzyme

organisms = get_organisms_for_enzyme("6.3.5.5")  # Glutamine synthetase
print(f"Found {len(organisms)} organisms with this enzyme")
```

## 5. Environmental Parameters

**Get pH and Temperature Data**:
```python
from scripts.brenda_queries import get_environmental_parameters

params = get_environmental_parameters("1.1.1.1")

print(f"Optimal pH range: {params['ph_range']}")
print(f"Optimal temperature: {params['optimal_temperature']}")
print(f"Stability pH: {params['stability_ph']}")
print(f"Temperature stability: {params['temperature_stability']}")
```

**Cofactor Requirements**:
```python
from scripts.brenda_queries import get_cofactor_requirements

cofactors = get_cofactor_requirements("1.1.1.1")
for cofactor in cofactors:
    print(f"Cofactor: {cofactor['name']}")
    print(f"Type: {cofactor['type']}")
    print(f"Concentration: {cofactor['concentration']}")
```

## 6. Substrate Specificity

**Get Substrate Specificity Data**:
```python
from scripts.brenda_queries import get_substrate_specificity

specificity = get_substrate_specificity("1.1.1.1")

for substrate in specificity:
    print(f"Substrate: {substrate['name']}")
    print(f"Km: {substrate['km']}")
    print(f"Vmax: {substrate['vmax']}")
    print(f"kcat: {substrate['kcat']}")
    print(f"Specificity constant: {substrate['kcat_km_ratio']}")
```

**Compare Substrate Preferences**:
```python
from scripts.brenda_queries import compare_substrate_affinity

comparison = compare_substrate_affinity("1.1.1.1")
sorted_by_km = sorted(comparison, key=lambda x: x['km'])

for substrate in sorted_by_km[:5]:  # Top 5 lowest Km
    print(f"{substrate['name']}: Km = {substrate['km']}")
```

## 7. Inhibition and Activation

**Get Inhibitor Information**:
```python
from scripts.brenda_queries import get_inhibitors

inhibitors = get_inhibitors("1.1.1.1")

for inhibitor in inhibitors:
    print(f"Inhibitor: {inhibitor['name']}")
    print(f"Type: {inhibitor['type']}")
    print(f"Ki: {inhibitor['ki']}")
    print(f"IC50: {inhibitor['ic50']}")
```

**Get Activator Information**:
```python
from scripts.brenda_queries import get_activators

activators = get_activators("1.1.1.1")

for activator in activators:
    print(f"Activator: {activator['name']}")
    print(f"Effect: {activator['effect']}")
    print(f"Mechanism: {activator['mechanism']}")
```

## 8. Enzyme Engineering Support

**Find Thermophilic Homologs**:
```python
from scripts.brenda_queries import find_thermophilic_homologs

thermophilic = find_thermophilic_homologs("1.1.1.1", min_temp=50)

for enzyme in thermophilic:
    print(f"Organism: {enzyme['organism']}")
    print(f"Optimal temp: {enzyme['optimal_temperature']}")
    print(f"Km: {enzyme['km']}")
```

**Find Alkaline/Acid Stable Variants**:
```python
from scripts.brenda_queries import find_ph_stable_variants

alkaline = find_ph_stable_variants("1.1.1.1", min_ph=8.0)
acidic = find_ph_stable_variants("1.1.1.1", max_ph=6.0)
```

## 9. Kinetic Modeling

**Get Kinetic Parameters for Modeling**:
```python
from scripts.brenda_queries import get_modeling_parameters

model_data = get_modeling_parameters("1.1.1.1", substrate="ethanol")

print(f"Km: {model_data['km']}")
print(f"Vmax: {model_data['vmax']}")
print(f"kcat: {model_data['kcat']}")
print(f"Enzyme concentration: {model_data['enzyme_conc']}")
print(f"Temperature: {model_data['temperature']}")
print(f"pH: {model_data['ph']}")
```

**Generate Michaelis-Menten Plots**:
```python
from scripts.brenda_visualization import plot_michaelis_menten

# Generate kinetic plots
plot_michaelis_menten("1.1.1.1", substrate="ethanol")
```

# BRENDA — Data Formats, Parsing & Best Practices

## BRENDA Response Format

BRENDA returns data in specific delimited formats that need parsing.

**Km Value Format**:
```
organism*Escherichia coli#substrate*ethanol#kmValue*1.2#kmValueMaximum*#commentary*pH 7.4, 25°C#ligandStructureId*#literature*
```

**Reaction Format**:
```
ecNumber*1.1.1.1#organism*Saccharomyces cerevisiae#reaction*ethanol + NAD+ <=> acetaldehyde + NADH + H+#commentary*#literature*
```

## Data Extraction Patterns

```python
import re

def parse_brenda_field(data, field_name):
    """Extract specific field from BRENDA data entry"""
    pattern = f"{field_name}\\*([^#]*)"
    match = re.search(pattern, data)
    return match.group(1) if match else None

def extract_multiple_values(data, field_name):
    """Extract multiple values for a field"""
    pattern = f"{field_name}\\*([^#]*)"
    matches = re.findall(pattern, data)
    return [match for match in matches if match.strip()]
```

## API Rate Limits and Best Practices

**Rate Limits**:
- BRENDA asks for no more than one request per second (the client enforces it)
- Faster clients may be identified as bots and blocked

**Best Practices**:
1. **Cache results**: Store frequently accessed enzyme data locally
2. **Batch queries**: Combine related requests when possible
3. **Use specific searches**: Narrow down by organism, substrate when possible
4. **Handle missing data**: Not all enzymes have complete data
5. **Validate EC numbers**: Ensure EC numbers are in correct format
6. **Implement delays**: Add delays between consecutive requests
7. **Use wildcards wisely**: Use '*' for broader searches when appropriate
8. **Monitor quota**: Track your API usage

**Error Handling**:
```python
from scripts.brenda_client import get_km_values, get_reactions
from zeep.exceptions import Fault, TransportError

try:
    km_data = get_km_values("1.1.1.1")
except RuntimeError as e:
    print(f"Authentication error: {e}")
except Fault as e:
    print(f"BRENDA API error: {e}")
except TransportError as e:
    print(f"Network error: {e}")
except Exception as e:
    print(f"Unexpected error: {e}")
```

## Troubleshooting

**Authentication Errors**:
- Verify BRENDA_EMAIL and BRENDA_PASSWORD in .env file
- Check for correct spelling (note BRENDA_EMIAL legacy support)
- Ensure BRENDA account is active and has API access

**No Results Returned**:
- Try broader searches with wildcards (*)
- Check EC number format (e.g., "1.1.1.1" not "1.1.1")
- Verify substrate spelling and naming
- Some enzymes may have limited data in BRENDA

**Rate Limiting**:
- Add delays between requests (0.5-1 second)
- Cache results locally
- Use more specific queries to reduce data volume
- Consider batch operations for multiple queries

**Network Errors**:
- Check internet connection
- BRENDA server may be temporarily unavailable
- Try again after a few minutes
- Consider using VPN if geo-restricted

**Data Format Issues**:
- Use the provided parsing functions in scripts
- BRENDA data can be inconsistent in formatting
- Handle missing fields gracefully
- Validate parsed data before use

**Performance Issues**:
- Large queries can be slow; limit search scope
- Use specific organism or substrate filters
- Consider asynchronous processing for batch operations
- Monitor memory usage with large datasets

## Additional Resources

- BRENDA Home: https://www.brenda-enzymes.org/
- BRENDA SOAP API Documentation: https://www.brenda-enzymes.org/soap.php
- Enzyme Commission (EC) Numbers: https://www.qmul.ac.uk/sbcs/iubmb/enzyme/
- Zeep SOAP Client: https://python-zeep.readthedocs.io/
- Enzyme Nomenclature: https://www.iubmb.org/enzyme/

# BRENDA — Helper Scripts Reference

This skill ships three helper scripts under `scripts/`. Function inventories below.

## scripts/brenda_queries.py

High-level functions for enzyme data analysis:

- `parse_km_entry(entry)`: Parse BRENDA Km data entries
- `parse_reaction_entry(entry)`: Parse reaction data entries
- `extract_organism_data(entry)`: Extract organism-specific information
- `search_enzymes_by_substrate(substrate, limit)`: Find enzymes for substrates
- `search_enzymes_by_product(product, limit)`: Find enzymes producing products
- `compare_across_organisms(ec_number, organisms)`: Compare enzyme properties
- `get_environmental_parameters(ec_number)`: Get pH and temperature data
- `get_cofactor_requirements(ec_number)`: Get cofactor information
- `get_substrate_specificity(ec_number)`: Analyze substrate preferences
- `get_inhibitors(ec_number)`: Get enzyme inhibition data
- `get_activators(ec_number)`: Get enzyme activation data
- `find_thermophilic_homologs(ec_number, min_temp)`: Find heat-stable variants
- `get_modeling_parameters(ec_number, substrate)`: Get parameters for kinetic modeling
- `export_kinetic_data(ec_number, format, filename)`: Export data to file

```python
from scripts.brenda_queries import search_enzymes_by_substrate, compare_across_organisms

# Search for enzymes
enzymes = search_enzymes_by_substrate("glucose", limit=20)

# Compare across organisms
comparison = compare_across_organisms("1.1.1.1", ["E. coli", "S. cerevisiae"])
```

## scripts/brenda_visualization.py

Visualization functions for enzyme data:

- `plot_kinetic_parameters(ec_number)`: Plot Km and kcat distributions
- `plot_organism_comparison(ec_number, organisms)`: Compare organisms
- `plot_pH_profiles(ec_number)`: Plot pH activity profiles
- `plot_temperature_profiles(ec_number)`: Plot temperature activity profiles
- `plot_substrate_specificity(ec_number)`: Visualize substrate preferences
- `plot_michaelis_menten(ec_number, substrate)`: Generate kinetic curves
- `create_heatmap_data(enzymes, parameters)`: Create data for heatmaps
- `generate_summary_plots(ec_number)`: Create comprehensive enzyme overview

```python
from scripts.brenda_visualization import plot_kinetic_parameters, plot_michaelis_menten

# Plot kinetic parameters
plot_kinetic_parameters("1.1.1.1")

# Generate Michaelis-Menten curve
plot_michaelis_menten("1.1.1.1", substrate="ethanol")
```

## scripts/enzyme_pathway_builder.py

Build enzymatic pathways and retrosynthetic routes:

- `find_pathway_for_product(product, max_steps)`: Find enzymatic pathways
- `build_retrosynthetic_tree(target, depth)`: Build retrosynthetic tree
- `suggest_enzyme_substitutions(ec_number, criteria)`: Suggest enzyme alternatives
- `calculate_pathway_feasibility(pathway)`: Evaluate pathway viability
- `optimize_pathway_conditions(pathway)`: Suggest optimal conditions
- `generate_pathway_report(pathway, filename)`: Create detailed pathway report

```python
from scripts.enzyme_pathway_builder import find_pathway_for_product, build_retrosynthetic_tree

# Find pathway to product
pathway = find_pathway_for_product("lactate", max_steps=3)

# Build retrosynthetic tree
tree = build_retrosynthetic_tree("lactate", depth=2)
```

# BRENDA — Common Workflows

End-to-end recipes built on the helper scripts. See `helper_scripts.md` for the function
inventory and `capabilities.md` for per-capability snippets.

## Workflow 1: Enzyme Discovery for New Substrate

Find suitable enzymes for a specific substrate:

Note: `getKmValue` entries carry no `ecNumber`, so to recover EC numbers and
reaction equations for a substrate, discover via reaction data (`search_by_pattern`),
then pull kinetics with `get_km_values`.

```python
from scripts.brenda_client import get_km_values
from scripts.brenda_queries import search_by_pattern

# Discover enzymes whose reactions mention the substrate (these entries have ecNumber)
substrate = "2-phenylethanol"
hits = search_by_pattern(substrate, limit=15)

print(f"Found {len(hits)} reaction hits for {substrate}")
for hit in hits:
    print(f"EC {hit['ec_number']} ({hit['organism']}): {hit['reaction']}")

# Get kinetic data for the first candidate with an EC number
candidates = [h for h in hits if h['ec_number']]
if candidates:
    best_ec = candidates[0]['ec_number']
    km_data = get_km_values(best_ec, substrate=substrate)

    if km_data:
        print(f"Kinetic data for {best_ec}:")
        for entry in km_data[:3]:  # First 3 entries
            print(f"  {entry}")
```

## Workflow 2: Cross-Organism Enzyme Comparison

Compare enzyme properties across different organisms:

```python
from scripts.brenda_queries import compare_across_organisms, get_environmental_parameters

# Define organisms for comparison
organisms = [
    "Escherichia coli",
    "Saccharomyces cerevisiae",
    "Bacillus subtilis",
    "Thermus thermophilus"
]

# Compare alcohol dehydrogenase
comparison = compare_across_organisms("1.1.1.1", organisms)

print("Cross-organism comparison:")
for org_data in comparison:
    print(f"\n{org_data['organism']}:")
    print(f"  Average Km: {org_data['average_km']}")
    print(f"  Optimal pH: {org_data['optimal_ph']}")
    print(f"  Temperature: {org_data['optimal_temperature']}°C")

# Get detailed environmental parameters
env_params = get_environmental_parameters("1.1.1.1")
print(f"\nOverall optimal pH range: {env_params['ph_range']}")
```

## Workflow 3: Enzyme Engineering Target Identification

Find engineering opportunities for enzyme improvement:

```python
from scripts.brenda_queries import (
    find_thermophilic_homologs,
    find_ph_stable_variants,
    compare_substrate_affinity
)

# Find thermophilic variants for heat stability
thermophilic = find_thermophilic_homologs("1.1.1.1", min_temp=50)
print(f"Found {len(thermophilic)} thermophilic variants")

# Find alkaline-stable variants
alkaline = find_ph_stable_variants("1.1.1.1", min_ph=8.0)
print(f"Found {len(alkaline)} alkaline-stable variants")

# Compare substrate specificities for engineering targets
specificity = compare_substrate_affinity("1.1.1.1")
print("Substrate affinity ranking:")
for i, sub in enumerate(specificity[:5]):
    print(f"  {i+1}. {sub['name']}: Km = {sub['km']}")
```

## Workflow 4: Enzymatic Pathway Construction

Build enzymatic synthesis pathways:

```python
from scripts.enzyme_pathway_builder import (
    find_pathway_for_product,
    build_retrosynthetic_tree,
    calculate_pathway_feasibility
)

# Find pathway to target product
target = "lactate"
pathway = find_pathway_for_product(target, max_steps=3)

if pathway:
    print(f"Found pathway to {target}:")
    for i, step in enumerate(pathway['steps']):
        print(f"  Step {i+1}: {step['reaction']}")
        print(f"    Enzyme: EC {step['ec_number']}")
        print(f"    Organism: {step['organism']}")

# Evaluate pathway feasibility
feasibility = calculate_pathway_feasibility(pathway)
print(f"\nPathway feasibility score: {feasibility['score']}/10")
print(f"Potential issues: {feasibility['warnings']}")
```

## Workflow 5: Kinetic Parameter Analysis

Comprehensive kinetic analysis for enzyme selection:

```python
from scripts.brenda_client import get_km_values
from scripts.brenda_queries import parse_km_entry, get_modeling_parameters
from scripts.brenda_visualization import plot_kinetic_parameters

# Get comprehensive kinetic data
ec_number = "1.1.1.1"
km_data = get_km_values(ec_number)

# Analyze kinetic parameters. parse_km_entry only sets 'km_value_numeric'
# when the raw kmValue field contained a parseable number, so filter on it.
all_entries = []
for entry in km_data:
    parsed = parse_km_entry(entry)
    if parsed.get('km_value_numeric') is not None:
        all_entries.append(parsed)

print(f"Analyzed {len(all_entries)} kinetic entries")

# Find best kinetic performer (lowest Km = highest affinity)
best_km = min(all_entries, key=lambda x: x['km_value_numeric'])
print(f"\nBest kinetic performer:")
print(f"  Organism: {best_km.get('organism')}")
print(f"  Substrate: {best_km.get('substrate')}")
print(f"  Km: {best_km['km_value_numeric']}")

# Get modeling parameters
model_data = get_modeling_parameters(ec_number, substrate=best_km['substrate'])
print(f"\nModeling parameters:")
print(f"  Km: {model_data['km']}")
print(f"  kcat: {model_data['kcat']}")
print(f"  Vmax: {model_data['vmax']}")

# Generate visualization
plot_kinetic_parameters(ec_number)
```

## Workflow 6: Industrial Enzyme Selection

Select enzymes for industrial applications:

```python
from scripts.brenda_queries import (
    find_thermophilic_homologs,
    get_environmental_parameters,
    get_inhibitors
)

# Industrial criteria: high temperature tolerance, organic solvent resistance
target_enzyme = "1.1.1.1"

# Find thermophilic variants
thermophilic = find_thermophilic_homologs(target_enzyme, min_temp=60)
print(f"Thermophilic candidates: {len(thermophilic)}")

# Check solvent tolerance (inhibitor data)
inhibitors = get_inhibitors(target_enzyme)
solvent_tolerant = [
    inv for inv in inhibitors
    if 'ethanol' not in inv['name'].lower() and
       'methanol' not in inv['name'].lower()
]

print(f"Solvent tolerant candidates: {len(solvent_tolerant)}")

# Evaluate top candidates
for candidate in thermophilic[:3]:
    print(f"\nCandidate: {candidate['organism']}")
    print(f"  Optimal temp: {candidate['optimal_temperature']}°C")
    print(f"  Km: {candidate['km']}")
    print(f"  pH range: {candidate.get('ph_range', 'N/A')}")
```

"""
BRENDA SOAP API client.

Thin wrapper around the official BRENDA SOAP API (https://www.brenda-enzymes.org/)
using the zeep SOAP client. Implements credential loading, SHA-256 password
hashing, a singleton SOAP client, and the convenience query functions used
throughout this skill.

Authentication:
    BRENDA requires a registered account. Provide credentials via a .env file
    (BRENDA_EMAIL / BRENDA_PASSWORD) or environment variables. The password is
    SHA-256 hashed before being sent, per the BRENDA SOAP specification.

Calling convention (brenda_zeep.wsdl):
    Every operation takes SEPARATE string arguments in WSDL order —
    email, password-hash, then "field*value" tokens — i.e.
    ``client.service.getKmValue(email, pw_hash, "ecNumber*1.1.1.1", ...)``.
    Passing one comma-joined string (the pre-zeep SOAPpy style) fails in zeep
    with "Missing element password". Results come back as lists of typed
    objects; split_entries() renders each one as the legacy
    "field*value#field*value" string that the parsers in brenda_queries expect.

Usage policy:
    BRENDA asks for at most one request per second; call_brenda() enforces it.
    Data are licensed CC BY 4.0 — cite BRENDA when you use them.

Installation:
    uv pip install zeep requests

Usage:
    from scripts.brenda_client import get_km_values, get_reactions

    km_data = get_km_values("1.1.1.1", organism="Saccharomyces cerevisiae")
    reactions = get_reactions("1.1.1.1")
"""

import hashlib
import os
import time
from pathlib import Path
from typing import List

from zeep import Client, Settings
from zeep.transports import Transport

WSDL_URL = "https://www.brenda-enzymes.org/soap/brenda_zeep.wsdl"

_CLIENT = None  # singleton zeep Client
_MIN_INTERVAL = 1.0  # seconds between calls (BRENDA: max one request per second)
_last_call = 0.0


def load_env_from_file(path: str = ".env") -> None:
    """Load KEY=VALUE pairs from a .env file into os.environ.

    Existing environment variables are not overwritten. Lines that are blank
    or start with '#' are ignored.
    """
    env_path = Path(path)
    if not env_path.exists():
        return
    for line in env_path.read_text().splitlines():
        line = line.strip()
        if not line or line.startswith("#") or "=" not in line:
            continue
        key, _, value = line.partition("=")
        key = key.strip()
        value = value.strip().strip('"').strip("'")
        os.environ.setdefault(key, value)


def _get_credentials() -> tuple:
    """Return (email, password) from the environment.

    Falls back to the legacy misspelled BRENDA_EMIAL variable for the email.
    Raises RuntimeError if either credential is missing.
    """
    load_env_from_file()
    email = os.environ.get("BRENDA_EMAIL") or os.environ.get("BRENDA_EMIAL")
    password = os.environ.get("BRENDA_PASSWORD")
    if not email or not password:
        raise RuntimeError(
            "BRENDA credentials missing. Set BRENDA_EMAIL and BRENDA_PASSWORD "
            "in your environment or a .env file."
        )
    return email, password


def _hash_password(password: str) -> str:
    """Return the SHA-256 hex digest of a plaintext password."""
    return hashlib.sha256(password.encode("utf-8")).hexdigest()


def _get_client() -> Client:
    """Initialize (once) and return the singleton zeep SOAP client."""
    global _CLIENT
    if _CLIENT is None:
        settings = Settings(strict=False, xml_huge_tree=True)
        transport = Transport(timeout=60)
        _CLIENT = Client(WSDL_URL, settings=settings, transport=transport)
    return _CLIENT


def call_brenda(action: str, parameters: List[str]):
    """Execute a BRENDA SOAP action.

    Args:
        action: SOAP method name, e.g. "getKmValue" or "getReaction".
        parameters: Field tokens in the operation's WSDL order, such as
            ["ecNumber*1.1.1.1", "organism*Homo sapiens", "kmValue*", ...].
            Email and the SHA-256-hashed password are prepended automatically.

    Returns:
        The raw zeep response (usually a list of typed result objects).
    """
    global _last_call
    email, password = _get_credentials()
    hashed = _hash_password(password)
    client = _get_client()
    wait = _MIN_INTERVAL - (time.monotonic() - _last_call)
    if wait > 0:
        time.sleep(wait)
    method = getattr(client.service, action)
    try:
        # zeep maps positional arguments onto the WSDL message parts in order.
        return method(email, hashed, *parameters)
    finally:
        _last_call = time.monotonic()


def _entry_to_string(item) -> str:
    """Render one zeep result object as a legacy 'field*value#field*value' string."""
    from zeep.helpers import serialize_object

    data = serialize_object(item)
    if not isinstance(data, dict):
        return str(data)
    parts = []
    for key, value in data.items():
        if value is None:
            value = ""
        elif isinstance(value, (list, tuple)):
            value = ", ".join(str(v) for v in value)
        parts.append(f"{key}*{value}")
    return "#".join(parts)


def split_entries(return_text) -> List[str]:
    """Normalize a BRENDA response into a list of entry strings.

    The zeep WSDL returns a list of typed objects; each is rendered as the
    legacy "field*value#field*value" string. Plain-string responses (older
    SOAP clients) separate records with '!'. Returns [] for empty input.
    """
    if not return_text:
        return []
    if isinstance(return_text, (list, tuple)):
        entries = [_entry_to_string(item) for item in return_text]
        return [entry for entry in entries if entry.strip()]
    return [entry for entry in str(return_text).split("!") if entry.strip()]


def get_km_values(ec_number: str, organism: str = "*", substrate: str = "*") -> List[str]:
    """Retrieve Km values for an enzyme.

    Args:
        ec_number: Enzyme Commission number (e.g., "1.1.1.1").
        organism: Organism name; "*" matches all organisms.
        substrate: Substrate name; "*" matches all substrates.

    Returns:
        List of raw BRENDA Km data entries.
    """
    # WSDL order: ecNumber, organism, kmValue, kmValueMaximum, substrate,
    # commentary, ligandStructureId, literature.
    parameters = [
        f"ecNumber*{ec_number}",
        f"organism*{'' if organism == '*' else organism}",
        "kmValue*",
        "kmValueMaximum*",
        f"substrate*{'' if substrate == '*' else substrate}",
        "commentary*",
        "ligandStructureId*",
        "literature*",
    ]
    return split_entries(call_brenda("getKmValue", parameters))


def get_reactions(ec_number: str, organism: str = "*", reaction: str = "*") -> List[str]:
    """Retrieve reaction data for an enzyme.

    Args:
        ec_number: Enzyme Commission number (e.g., "1.1.1.1").
        organism: Organism name; "*" matches all organisms.
        reaction: Reaction pattern; "*" matches all reactions.

    Returns:
        List of raw BRENDA reaction data entries.
    """
    parameters = [
        f"ecNumber*{ec_number}",
        f"organism*{'' if organism == '*' else organism}",
        f"reaction*{'' if reaction == '*' else reaction}",
        "commentary*",
        "literature*",
    ]
    return split_entries(call_brenda("getReaction", parameters))

"""
BRENDA Database Query Utilities

This module provides high-level functions for querying and analyzing
enzyme data from the BRENDA database using the SOAP API.

Key features:
- Parse BRENDA response data entries
- Search for enzymes by substrate/product
- Compare enzyme properties across organisms
- Retrieve kinetic parameters and environmental conditions
- Analyze substrate specificity and inhibition
- Support for enzyme engineering and pathway design
- Export data in various formats

Installation:
    uv pip install zeep requests pandas

Usage:
    from scripts.brenda_queries import search_enzymes_by_substrate, compare_across_organisms

    enzymes = search_enzymes_by_substrate("glucose", limit=20)
    comparison = compare_across_organisms("1.1.1.1", ["E. coli", "S. cerevisiae"])
"""

import re
import time
import json
from typing import List, Dict, Any

try:
    from zeep import Client, Settings  # noqa: F401  # availability probe for ZEEP_AVAILABLE
    from zeep.exceptions import Fault, TransportError  # noqa: F401  # availability probe
    ZEEP_AVAILABLE = True
except ImportError:
    print("Warning: zeep not installed. Install with: uv pip install zeep")
    ZEEP_AVAILABLE = False

try:
    import requests  # noqa: F401  # availability probe for REQUESTS_AVAILABLE
    REQUESTS_AVAILABLE = True
except ImportError:
    print("Warning: requests not installed. Install with: uv pip install requests")
    REQUESTS_AVAILABLE = False

try:
    import pandas as pd
    PANDAS_AVAILABLE = True
except ImportError:
    print("Warning: pandas not installed. Install with: uv pip install pandas")
    PANDAS_AVAILABLE = False

# Import the brenda_client shipped alongside this module.
try:
    from scripts.brenda_client import get_km_values, get_reactions, call_brenda
    BRENDA_CLIENT_AVAILABLE = True
except ImportError:
    try:
        from brenda_client import get_km_values, get_reactions, call_brenda  # noqa: F401  # module-path fallback for the package-path import above
        BRENDA_CLIENT_AVAILABLE = True
    except ImportError:
        print("Warning: brenda_client not available")
        BRENDA_CLIENT_AVAILABLE = False


def validate_dependencies():
    """Validate that required dependencies are installed."""
    missing = []
    if not ZEEP_AVAILABLE:
        missing.append("zeep")
    if not REQUESTS_AVAILABLE:
        missing.append("requests")
    if not BRENDA_CLIENT_AVAILABLE:
        missing.append("brenda_client")
    if missing:
        raise ImportError(f"Missing required dependencies: {', '.join(missing)}")


def parse_km_entry(entry: str) -> Dict[str, Any]:
    """Parse a BRENDA Km value entry into structured data."""
    if not entry or not isinstance(entry, str):
        return {}

    parsed = {}
    parts = entry.split('#')

    for part in parts:
        if '*' in part:
            key, value = part.split('*', 1)
            parsed[key.strip()] = value.strip()

    # Extract numeric values from kmValue
    if 'kmValue' in parsed:
        km_value = parsed['kmValue']
        # Extract first numeric value (in mM typically)
        numeric_match = re.search(r'(\d+\.?\d*)', km_value)
        if numeric_match:
            parsed['km_value_numeric'] = float(numeric_match.group(1))

    # Extract pH from commentary
    if 'commentary' in parsed:
        commentary = parsed['commentary']
        ph_match = re.search(r'pH\s*([0-9.]+)', commentary)
        if ph_match:
            parsed['ph'] = float(ph_match.group(1))

        temp_match = re.search(r'(\d+)\s*°?C', commentary)
        if temp_match:
            parsed['temperature'] = float(temp_match.group(1))

    return parsed


def parse_reaction_entry(entry: str) -> Dict[str, Any]:
    """Parse a BRENDA reaction entry into structured data."""
    if not entry or not isinstance(entry, str):
        return {}

    parsed = {}
    parts = entry.split('#')

    for part in parts:
        if '*' in part:
            key, value = part.split('*', 1)
            parsed[key.strip()] = value.strip()

    # Parse reaction equation
    if 'reaction' in parsed:
        reaction = parsed['reaction']
        # Extract reactants and products
        if '<=>' in reaction:
            reactants, products = reaction.split('<=>', 1)
        elif '->' in reaction:
            reactants, products = reaction.split('->', 1)
        elif '=' in reaction:
            reactants, products = reaction.split('=', 1)
        else:
            reactants, products = reaction, ''

        parsed['reactants'] = [r.strip() for r in reactants.split('+')]
        parsed['products'] = [p.strip() for p in products.split('+')]

    return parsed


def extract_organism_data(entry: str) -> Dict[str, Any]:
    """Extract organism-specific information from BRENDA entry."""
    parsed = parse_km_entry(entry) if 'kmValue' in entry else parse_reaction_entry(entry)

    if 'organism' in parsed:
        return {
            'organism': parsed['organism'],
            'ec_number': parsed.get('ecNumber', ''),
            'substrate': parsed.get('substrate', ''),
            'km_value': parsed.get('kmValue', ''),
            'km_numeric': parsed.get('km_value_numeric', None),
            'ph': parsed.get('ph', None),
            'temperature': parsed.get('temperature', None),
            'commentary': parsed.get('commentary', ''),
            'literature': parsed.get('literature', '')
        }

    return {}


def search_enzymes_by_substrate(substrate: str, limit: int = 50) -> List[Dict[str, Any]]:
    """Search for enzymes that act on a specific substrate."""
    validate_dependencies()

    enzymes = []

    # Search for Km values with the substrate
    try:
        km_data = get_km_values("*", substrate=substrate)
        time.sleep(0.5)  # Rate limiting

        for entry in km_data[:limit]:
            parsed = parse_km_entry(entry)
            if parsed:
                enzymes.append({
                    'ec_number': parsed.get('ecNumber', ''),
                    'organism': parsed.get('organism', ''),
                    'substrate': parsed.get('substrate', ''),
                    'km_value': parsed.get('kmValue', ''),
                    'km_numeric': parsed.get('km_value_numeric', None),
                    'commentary': parsed.get('commentary', '')
                })
    except Exception as e:
        print(f"Error searching enzymes by substrate: {e}")

    # Remove duplicates based on EC number and organism
    unique_enzymes = []
    seen = set()
    for enzyme in enzymes:
        key = (enzyme['ec_number'], enzyme['organism'])
        if key not in seen:
            seen.add(key)
            unique_enzymes.append(enzyme)

    return unique_enzymes[:limit]


def search_enzymes_by_product(product: str, limit: int = 50) -> List[Dict[str, Any]]:
    """Search for enzymes that produce a specific product."""
    validate_dependencies()

    enzymes = []

    # Search for reactions containing the product
    try:
        # This is a simplified approach - in practice you might need
        # more sophisticated pattern matching for products
        reactions = get_reactions("*", reaction=f"*{product}*")
        time.sleep(0.5)  # Rate limiting

        for entry in reactions[:limit]:
            parsed = parse_reaction_entry(entry)
            if parsed and 'products' in parsed:
                # Check if our target product is in the products list
                if any(product.lower() in prod.lower() for prod in parsed['products']):
                    enzymes.append({
                        'ec_number': parsed.get('ecNumber', ''),
                        'organism': parsed.get('organism', ''),
                        'reaction': parsed.get('reaction', ''),
                        'reactants': parsed.get('reactants', []),
                        'products': parsed.get('products', []),
                        'commentary': parsed.get('commentary', '')
                    })
    except Exception as e:
        print(f"Error searching enzymes by product: {e}")

    return enzymes[:limit]


def compare_across_organisms(ec_number: str, organisms: List[str]) -> List[Dict[str, Any]]:
    """Compare enzyme properties across different organisms."""
    validate_dependencies()

    comparison = []

    for organism in organisms:
        try:
            # Get Km data for this organism
            km_data = get_km_values(ec_number, organism=organism)
            time.sleep(0.5)  # Rate limiting

            if km_data:
                # Calculate statistics
                numeric_kms = []
                phs = []
                temperatures = []

                for entry in km_data:
                    parsed = parse_km_entry(entry)
                    if 'km_value_numeric' in parsed:
                        numeric_kms.append(parsed['km_value_numeric'])
                    if 'ph' in parsed:
                        phs.append(parsed['ph'])
                    if 'temperature' in parsed:
                        temperatures.append(parsed['temperature'])

                org_data = {
                    'organism': organism,
                    'ec_number': ec_number,
                    'data_points': len(km_data),
                    'average_km': sum(numeric_kms) / len(numeric_kms) if numeric_kms else None,
                    'min_km': min(numeric_kms) if numeric_kms else None,
                    'max_km': max(numeric_kms) if numeric_kms else None,
                    'optimal_ph': sum(phs) / len(phs) if phs else None,
                    'optimal_temperature': sum(temperatures) / len(temperatures) if temperatures else None,
                    'temperature_range': (min(temperatures), max(temperatures)) if temperatures else None
                }

                comparison.append(org_data)
            else:
                comparison.append({
                    'organism': organism,
                    'ec_number': ec_number,
                    'data_points': 0,
                    'note': 'No data found'
                })

        except Exception as e:
            print(f"Error comparing organism {organism}: {e}")
            comparison.append({
                'organism': organism,
                'ec_number': ec_number,
                'error': str(e)
            })

    return comparison


def get_organisms_for_enzyme(ec_number: str) -> List[str]:
    """Get list of organisms that have data for a specific enzyme."""
    validate_dependencies()

    try:
        km_data = get_km_values(ec_number)
        time.sleep(0.5)  # Rate limiting

        organisms = set()
        for entry in km_data:
            parsed = parse_km_entry(entry)
            if 'organism' in parsed:
                organisms.add(parsed['organism'])

        return sorted(list(organisms))

    except Exception as e:
        print(f"Error getting organisms for enzyme {ec_number}: {e}")
        return []


def get_environmental_parameters(ec_number: str) -> Dict[str, Any]:
    """Get environmental parameters (pH, temperature) for an enzyme."""
    validate_dependencies()

    try:
        km_data = get_km_values(ec_number)
        time.sleep(0.5)  # Rate limiting

        phs = []
        temperatures = []
        ph_stabilities = []
        temp_stabilities = []

        for entry in km_data:
            parsed = parse_km_entry(entry)

            if 'ph' in parsed:
                phs.append(parsed['ph'])
            if 'temperature' in parsed:
                temperatures.append(parsed['temperature'])

            # Check commentary for stability information
            commentary = parsed.get('commentary', '').lower()
            if 'stable' in commentary and 'ph' in commentary:
                # Extract pH stability range
                ph_range_match = re.search(r'ph\s*([\d.]+)\s*[-–]\s*([\d.]+)', commentary)
                if ph_range_match:
                    ph_stabilities.append((float(ph_range_match.group(1)), float(ph_range_match.group(2))))

            if 'stable' in commentary and ('temp' in commentary or '°c' in commentary):
                # Extract temperature stability
                temp_match = re.search(r'(\d+)\s*[-–]\s*(\d+)\s*°?c', commentary)
                if temp_match:
                    temp_stabilities.append((int(temp_match.group(1)), int(temp_match.group(2))))

        params = {
            'ec_number': ec_number,
            'data_points': len(km_data),
            'ph_range': (min(phs), max(phs)) if phs else None,
            'optimal_ph': sum(phs) / len(phs) if phs else None,
            'optimal_temperature': sum(temperatures) / len(temperatures) if temperatures else None,
            'temperature_range': (min(temperatures), max(temperatures)) if temperatures else None,
            'stability_ph': ph_stabilities[0] if ph_stabilities else None,
            'temperature_stability': temp_stabilities[0] if temp_stabilities else None
        }

        return params

    except Exception as e:
        print(f"Error getting environmental parameters for {ec_number}: {e}")
        return {'ec_number': ec_number, 'error': str(e)}


def get_cofactor_requirements(ec_number: str) -> List[Dict[str, Any]]:
    """Get cofactor requirements for an enzyme from reaction data."""
    validate_dependencies()

    cofactors = []

    try:
        reactions = get_reactions(ec_number)
        time.sleep(0.5)  # Rate limiting

        for entry in reactions:
            parsed = parse_reaction_entry(entry)
            if parsed and 'reactants' in parsed:
                # Look for common cofactors in reactants
                common_cofactors = [
                    'NAD+', 'NADH', 'NADP+', 'NADPH',
                    'ATP', 'ADP', 'AMP',
                    'FAD', 'FADH2',
                    'CoA', 'acetyl-CoA',
                    'pyridoxal phosphate', 'PLP',
                    'biotin',
                    'heme', 'iron-sulfur'
                ]

                for reactant in parsed['reactants']:
                    for cofactor in common_cofactors:
                        if cofactor.lower() in reactant.lower():
                            cofactors.append({
                                'name': cofactor,
                                'full_name': reactant,
                                'type': 'oxidoreductase' if 'NAD' in cofactor else 'other',
                                'organism': parsed.get('organism', ''),
                                'ec_number': ec_number
                            })

    except Exception as e:
        print(f"Error getting cofactor requirements for {ec_number}: {e}")

    # Remove duplicates
    unique_cofactors = []
    seen = set()
    for cofactor in cofactors:
        key = (cofactor['name'], cofactor['organism'])
        if key not in seen:
            seen.add(key)
            unique_cofactors.append(cofactor)

    return unique_cofactors


def get_substrate_specificity(ec_number: str) -> List[Dict[str, Any]]:
    """Get substrate specificity data for an enzyme."""
    validate_dependencies()

    specificity = []

    try:
        km_data = get_km_values(ec_number)
        time.sleep(0.5)  # Rate limiting

        substrate_data = {}

        for entry in km_data:
            parsed = parse_km_entry(entry)
            if 'substrate' in parsed and 'km_value_numeric' in parsed:
                substrate = parsed['substrate']
                if substrate not in substrate_data:
                    substrate_data[substrate] = {
                        'name': substrate,
                        'km_values': [],
                        'organisms': set(),
                        'vmax_values': [],  # If available
                        'kcat_values': []   # If available
                    }

                substrate_data[substrate]['km_values'].append(parsed['km_value_numeric'])
                if 'organism' in parsed:
                    substrate_data[substrate]['organisms'].add(parsed['organism'])

        # Calculate summary statistics
        for substrate, data in substrate_data.items():
            if data['km_values']:
                specificity.append({
                    'name': substrate,
                    'km': sum(data['km_values']) / len(data['km_values']),
                    'min_km': min(data['km_values']),
                    'max_km': max(data['km_values']),
                    'data_points': len(data['km_values']),
                    'organisms': list(data['organisms']),
                    'vmax': sum(data['vmax_values']) / len(data['vmax_values']) if data['vmax_values'] else None,
                    'kcat': sum(data['kcat_values']) / len(data['kcat_values']) if data['kcat_values'] else None,
                    'kcat_km_ratio': None  # Would need kcat data to calculate
                })

        # Sort by Km (lower is better affinity)
        specificity.sort(key=lambda x: x['km'] if x['km'] else float('inf'))

    except Exception as e:
        print(f"Error getting substrate specificity for {ec_number}: {e}")

    return specificity


def compare_substrate_affinity(ec_number: str) -> List[Dict[str, Any]]:
    """Compare substrate affinity for an enzyme."""
    return get_substrate_specificity(ec_number)


def get_inhibitors(ec_number: str) -> List[Dict[str, Any]]:
    """Get inhibitor information for an enzyme (from commentary)."""
    validate_dependencies()

    inhibitors = []

    try:
        km_data = get_km_values(ec_number)
        time.sleep(0.5)  # Rate limiting

        for entry in km_data:
            parsed = parse_km_entry(entry)
            commentary = parsed.get('commentary', '').lower()

            # Look for inhibitor keywords
            inhibitor_keywords = ['inhibited', 'inhibition', 'blocked', 'prevented', 'reduced']
            if any(keyword in commentary for keyword in inhibitor_keywords):
                # Try to extract inhibitor names (this is approximate)
                # Common inhibitors
                common_inhibitors = [
                    'iodoacetate', 'n-ethylmaleimide', 'p-chloromercuribenzoate',
                    'heavy metals', 'mercury', 'copper', 'zinc',
                    'cyanide', 'azide', 'carbon monoxide',
                    'edta', 'egta'
                ]

                for inhibitor in common_inhibitors:
                    if inhibitor in commentary:
                        inhibitors.append({
                            'name': inhibitor,
                            'type': 'irreversible' if 'iodoacetate' in inhibitor or 'maleimide' in inhibitor else 'reversible',
                            'organism': parsed.get('organism', ''),
                            'ec_number': ec_number,
                            'commentary': parsed.get('commentary', '')
                        })

    except Exception as e:
        print(f"Error getting inhibitors for {ec_number}: {e}")

    # Remove duplicates
    unique_inhibitors = []
    seen = set()
    for inhibitor in inhibitors:
        key = (inhibitor['name'], inhibitor['organism'])
        if key not in seen:
            seen.add(key)
            unique_inhibitors.append(inhibitor)

    return unique_inhibitors


def get_activators(ec_number: str) -> List[Dict[str, Any]]:
    """Get activator information for an enzyme (from commentary)."""
    validate_dependencies()

    activators = []

    try:
        km_data = get_km_values(ec_number)
        time.sleep(0.5)  # Rate limiting

        for entry in km_data:
            parsed = parse_km_entry(entry)
            commentary = parsed.get('commentary', '').lower()

            # Look for activator keywords
            activator_keywords = ['activated', 'stimulated', 'enhanced', 'increased']
            if any(keyword in commentary for keyword in activator_keywords):
                # Try to extract activator names (this is approximate)
                common_activators = [
                    'mg2+', 'mn2+', 'ca2+', 'zn2+',
                    'k+', 'na+',
                    'phosphate', 'pyrophosphate',
                    'dithiothreitol', 'dtt',
                    'β-mercaptoethanol'
                ]

                for activator in common_activators:
                    if activator in commentary:
                        activators.append({
                            'name': activator,
                            'type': 'metal ion' if '+' in activator else 'reducing agent' if 'dtt' in activator.lower() or 'mercapto' in activator.lower() else 'other',
                            'mechanism': 'allosteric' if 'allosteric' in commentary else 'cofactor' if 'cofactor' in commentary else 'unknown',
                            'organism': parsed.get('organism', ''),
                            'ec_number': ec_number,
                            'commentary': parsed.get('commentary', '')
                        })

    except Exception as e:
        print(f"Error getting activators for {ec_number}: {e}")

    # Remove duplicates
    unique_activators = []
    seen = set()
    for activator in activators:
        key = (activator['name'], activator['organism'])
        if key not in seen:
            seen.add(key)
            unique_activators.append(activator)

    return unique_activators


def find_thermophilic_homologs(ec_number: str, min_temp: int = 50) -> List[Dict[str, Any]]:
    """Find thermophilic homologs of an enzyme."""
    validate_dependencies()

    thermophilic = []

    try:
        organisms = get_organisms_for_enzyme(ec_number)

        for organism in organisms:
            # Check if organism might be thermophilic based on name
            thermophilic_keywords = ['therm', 'hypertherm', 'pyro']
            if any(keyword in organism.lower() for keyword in thermophilic_keywords):
                # Get kinetic data to extract temperature information
                km_data = get_km_values(ec_number, organism=organism)
                time.sleep(0.2)  # Rate limiting

                temperatures = []
                kms = []

                for entry in km_data:
                    parsed = parse_km_entry(entry)
                    if 'temperature' in parsed:
                        temperatures.append(parsed['temperature'])
                    if 'km_value_numeric' in parsed:
                        kms.append(parsed['km_value_numeric'])

                if temperatures and max(temperatures) >= min_temp:
                    thermophilic.append({
                        'organism': organism,
                        'ec_number': ec_number,
                        'optimal_temperature': max(temperatures),
                        'temperature_range': (min(temperatures), max(temperatures)),
                        'km': sum(kms) / len(kms) if kms else None,
                        'data_points': len(km_data)
                    })

    except Exception as e:
        print(f"Error finding thermophilic homologs for {ec_number}: {e}")

    return thermophilic


def find_ph_stable_variants(ec_number: str, min_ph: float = 8.0, max_ph: float = 6.0) -> List[Dict[str, Any]]:
    """Find pH-stable variants of an enzyme."""
    validate_dependencies()

    ph_stable = []

    try:
        organisms = get_organisms_for_enzyme(ec_number)

        for organism in organisms:
            km_data = get_km_values(ec_number, organism=organism)
            time.sleep(0.2)  # Rate limiting

            phs = []
            kms = []

            for entry in km_data:
                parsed = parse_km_entry(entry)
                if 'ph' in parsed:
                    phs.append(parsed['ph'])
                if 'km_value_numeric' in parsed:
                    kms.append(parsed['km_value_numeric'])

            if phs:
                ph_range = (min(phs), max(phs))
                is_alkaline_stable = min_ph and ph_range[0] >= min_ph
                is_acid_stable = max_ph and ph_range[1] <= max_ph

                if is_alkaline_stable or is_acid_stable:
                    ph_stable.append({
                        'organism': organism,
                        'ec_number': ec_number,
                        'ph_range': ph_range,
                        'optimal_ph': sum(phs) / len(phs),
                        'km': sum(kms) / len(kms) if kms else None,
                        'stability_type': 'alkaline' if is_alkaline_stable else 'acidic',
                        'data_points': len(km_data)
                    })

    except Exception as e:
        print(f"Error finding pH-stable variants for {ec_number}: {e}")

    return ph_stable


def get_modeling_parameters(ec_number: str, substrate: str = None) -> Dict[str, Any]:
    """Get parameters suitable for kinetic modeling."""
    validate_dependencies()

    try:
        if substrate:
            km_data = get_km_values(ec_number, substrate=substrate)
        else:
            km_data = get_km_values(ec_number)

        time.sleep(0.5)  # Rate limiting

        if not km_data:
            return {'ec_number': ec_number, 'error': 'No kinetic data found'}

        # Extract modeling parameters
        kms = []
        phs = []
        temperatures = []
        v_max_values = []
        kcat_values = []

        for entry in km_data:
            parsed = parse_km_entry(entry)

            if 'km_value_numeric' in parsed:
                kms.append(parsed['km_value_numeric'])
            if 'ph' in parsed:
                phs.append(parsed['ph'])
            if 'temperature' in parsed:
                temperatures.append(parsed['temperature'])

            # Look for Vmax and kcat in commentary (rare in BRENDA)
            commentary = parsed.get('commentary', '').lower()
            vmax_match = re.search(r'vmax\s*=\s*([\d.]+)', commentary)
            if vmax_match:
                v_max_values.append(float(vmax_match.group(1)))

            kcat_match = re.search(r'kcat\s*=\s*([\d.]+)', commentary)
            if kcat_match:
                kcat_values.append(float(kcat_match.group(1)))

        modeling_data = {
            'ec_number': ec_number,
            'substrate': substrate if substrate else 'various',
            'km': sum(kms) / len(kms) if kms else None,
            'km_std': (sum((x - sum(kms)/len(kms))**2 for x in kms) / len(kms))**0.5 if kms else None,
            'vmax': sum(v_max_values) / len(v_max_values) if v_max_values else None,
            'kcat': sum(kcat_values) / len(kcat_values) if kcat_values else None,
            'optimal_ph': sum(phs) / len(phs) if phs else None,
            'optimal_temperature': sum(temperatures) / len(temperatures) if temperatures else None,
            'data_points': len(km_data),
            'temperature': sum(temperatures) / len(temperatures) if temperatures else 25.0,  # Default to 25°C
            'ph': sum(phs) / len(phs) if phs else 7.0,  # Default to pH 7.0
            'enzyme_conc': 1.0,  # Default enzyme concentration (μM)
            'substrate_conc': None,  # Would be set by user
        }

        return modeling_data

    except Exception as e:
        return {'ec_number': ec_number, 'error': str(e)}


def export_kinetic_data(ec_number: str, format: str = 'csv', filename: str = None) -> str:
    """Export kinetic data to file."""
    validate_dependencies()

    if not filename:
        filename = f"brenda_kinetic_data_{ec_number.replace('.', '_')}.{format}"

    try:
        # Get all kinetic data
        km_data = get_km_values(ec_number)
        time.sleep(0.5)  # Rate limiting

        if not km_data:
            print(f"No kinetic data found for EC {ec_number}")
            return filename

        # Parse all entries
        parsed_data = []
        for entry in km_data:
            parsed = parse_km_entry(entry)
            if parsed:
                parsed_data.append(parsed)

        # Export based on format
        if format.lower() == 'csv':
            if parsed_data:
                df = pd.DataFrame(parsed_data)
                df.to_csv(filename, index=False)
            else:
                with open(filename, 'w', newline='') as f:
                    f.write('No data found')

        elif format.lower() == 'json':
            with open(filename, 'w') as f:
                json.dump(parsed_data, f, indent=2, default=str)

        elif format.lower() == 'excel':
            if parsed_data and PANDAS_AVAILABLE:
                df = pd.DataFrame(parsed_data)
                df.to_excel(filename, index=False)
            else:
                print("pandas required for Excel export")
                return filename

        print(f"Exported {len(parsed_data)} entries to {filename}")
        return filename

    except Exception as e:
        print(f"Error exporting data: {e}")
        return filename


def search_by_pattern(pattern: str, limit: int = 50) -> List[Dict[str, Any]]:
    """Search enzymes using a reaction pattern or keyword."""
    validate_dependencies()

    enzymes = []

    try:
        # Search reactions containing the pattern
        reactions = get_reactions("*", reaction=f"*{pattern}*")
        time.sleep(0.5)  # Rate limiting

        for entry in reactions[:limit]:
            parsed = parse_reaction_entry(entry)
            if parsed:
                enzymes.append({
                    'ec_number': parsed.get('ecNumber', ''),
                    'organism': parsed.get('organism', ''),
                    'reaction': parsed.get('reaction', ''),
                    'reactants': parsed.get('reactants', []),
                    'products': parsed.get('products', []),
                    'commentary': parsed.get('commentary', '')
                })

    except Exception as e:
        print(f"Error searching by pattern '{pattern}': {e}")

    return enzymes


if __name__ == "__main__":
    # Example usage
    print("BRENDA Database Query Examples")
    print("=" * 40)

    try:
        # Example 1: Search enzymes by substrate
        print("\n1. Searching enzymes for 'glucose':")
        enzymes = search_enzymes_by_substrate("glucose", limit=5)
        for enzyme in enzymes:
            print(f"  EC {enzyme['ec_number']}: {enzyme['organism']}")
            print(f"    Km: {enzyme['km_value']}")

        # Example 2: Compare across organisms
        print("\n2. Comparing alcohol dehydrogenase (1.1.1.1) across organisms:")
        organisms = ["Escherichia coli", "Saccharomyces cerevisiae", "Homo sapiens"]
        comparison = compare_across_organisms("1.1.1.1", organisms)
        for comp in comparison:
            if comp.get('data_points', 0) > 0:
                print(f"  {comp['organism']}:")
                print(f"    Avg Km: {comp.get('average_km', 'N/A')}")
                print(f"    Optimal pH: {comp.get('optimal_ph', 'N/A')}")

        # Example 3: Get environmental parameters
        print("\n3. Environmental parameters for 1.1.1.1:")
        params = get_environmental_parameters("1.1.1.1")
        if params.get('data_points', 0) > 0:
            print(f"  pH range: {params.get('ph_range', 'N/A')}")
            print(f"  Temperature range: {params.get('temperature_range', 'N/A')}")

    except Exception as e:
        print(f"Example failed: {e}")

"""
BRENDA Database Visualization Utilities

This module provides visualization functions for BRENDA enzyme data,
including kinetic parameters, environmental conditions, and pathway analysis.

Key features:
- Plot Km, kcat, and Vmax distributions
- Compare enzyme properties across organisms
- Visualize pH and temperature activity profiles
- Plot substrate specificity and affinity data
- Generate Michaelis-Menten curves
- Create heatmaps and correlation plots
- Support for pathway visualization

Installation:
    uv pip install matplotlib seaborn pandas numpy

Usage:
    from scripts.brenda_visualization import plot_kinetic_parameters, plot_michaelis_menten

    plot_kinetic_parameters("1.1.1.1")
    plot_michaelis_menten("1.1.1.1", substrate="ethanol")
"""

import numpy as np
from typing import List, Dict, Any
import matplotlib.pyplot as plt
import seaborn as sns
from pathlib import Path

try:
    import pandas as pd
    PANDAS_AVAILABLE = True
except ImportError:
    print("Warning: pandas not installed. Install with: uv pip install pandas")
    PANDAS_AVAILABLE = False

try:
    from brenda_queries import (
        get_km_values, parse_km_entry,
        compare_across_organisms,
        get_substrate_specificity, get_modeling_parameters,
    )
    BRENDA_QUERIES_AVAILABLE = True
except ImportError:
    print("Warning: brenda_queries not available")
    BRENDA_QUERIES_AVAILABLE = False


# Set style for plots
plt.style.use('default')
sns.set_palette("husl")


def validate_dependencies():
    """Validate that required dependencies are installed."""
    missing = []
    if not PANDAS_AVAILABLE:
        missing.append("pandas")
    if not BRENDA_QUERIES_AVAILABLE:
        missing.append("brenda_queries")
    if missing:
        raise ImportError(f"Missing required dependencies: {', '.join(missing)}")


def plot_kinetic_parameters(ec_number: str, save_path: str = None, show_plot: bool = True) -> str:
    """Plot kinetic parameter distributions for an enzyme."""
    validate_dependencies()

    try:
        # Get Km data
        km_data = get_km_values(ec_number)

        if not km_data:
            print(f"No kinetic data found for EC {ec_number}")
            return save_path

        # Parse data
        parsed_entries = []
        for entry in km_data:
            parsed = parse_km_entry(entry)
            if 'km_value_numeric' in parsed:
                parsed_entries.append(parsed)

        if not parsed_entries:
            print(f"No numeric Km data found for EC {ec_number}")
            return save_path

        # Create figure with subplots
        fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
        fig.suptitle(f'Kinetic Parameters for EC {ec_number}', fontsize=16, fontweight='bold')

        # Extract data
        km_values = [entry['km_value_numeric'] for entry in parsed_entries]
        organisms = [entry.get('organism', 'Unknown') for entry in parsed_entries]
        substrates = [entry.get('substrate', 'Unknown') for entry in parsed_entries]

        # Plot 1: Km distribution histogram
        ax1.hist(km_values, bins=30, alpha=0.7, edgecolor='black')
        ax1.set_xlabel('Km (mM)')
        ax1.set_ylabel('Frequency')
        ax1.set_title('Km Value Distribution')
        ax1.axvline(np.mean(km_values), color='red', linestyle='--', label=f'Mean: {np.mean(km_values):.2f}')
        ax1.axvline(np.median(km_values), color='blue', linestyle='--', label=f'Median: {np.median(km_values):.2f}')
        ax1.legend()

        # Plot 2: Km by organism (top 10)
        if PANDAS_AVAILABLE:
            df = pd.DataFrame({'Km': km_values, 'Organism': organisms})
            organism_means = df.groupby('Organism')['Km'].mean().sort_values(ascending=False).head(10)

            organism_means.plot(kind='bar', ax=ax2)
            ax2.set_ylabel('Mean Km (mM)')
            ax2.set_title('Mean Km by Organism (Top 10)')
            ax2.tick_params(axis='x', rotation=45)

        # Plot 3: Km by substrate (top 10)
        if PANDAS_AVAILABLE:
            df = pd.DataFrame({'Km': km_values, 'Substrate': substrates})
            substrate_means = df.groupby('Substrate')['Km'].mean().sort_values(ascending=False).head(10)

            substrate_means.plot(kind='bar', ax=ax3)
            ax3.set_ylabel('Mean Km (mM)')
            ax3.set_title('Mean Km by Substrate (Top 10)')
            ax3.tick_params(axis='x', rotation=45)

        # Plot 4: Box plot by organism (top 5)
        if PANDAS_AVAILABLE:
            top_organisms = df.groupby('Organism')['Km'].count().sort_values(ascending=False).head(5).index
            top_data = df[df['Organism'].isin(top_organisms)]

            sns.boxplot(data=top_data, x='Organism', y='Km', ax=ax4)
            ax4.set_ylabel('Km (mM)')
            ax4.set_title('Km Distribution by Organism (Top 5)')
            ax4.tick_params(axis='x', rotation=45)

        plt.tight_layout()

        # Save plot
        if save_path:
            plt.savefig(save_path, dpi=300, bbox_inches='tight')
            print(f"Kinetic parameters plot saved to {save_path}")

        if show_plot:
            plt.show()
        else:
            plt.close()

        return save_path or f"kinetic_parameters_{ec_number.replace('.', '_')}.png"

    except Exception as e:
        print(f"Error plotting kinetic parameters: {e}")
        return save_path


def plot_organism_comparison(ec_number: str, organisms: List[str], save_path: str = None, show_plot: bool = True) -> str:
    """Compare enzyme properties across multiple organisms."""
    validate_dependencies()

    try:
        # Get comparison data
        comparison = compare_across_organisms(ec_number, organisms)

        if not comparison:
            print(f"No comparison data found for EC {ec_number}")
            return save_path

        # Filter out entries with no data
        valid_data = [c for c in comparison if c.get('data_points', 0) > 0]

        if not valid_data:
            print(f"No valid data for organism comparison of EC {ec_number}")
            return save_path

        # Create figure
        fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
        fig.suptitle(f'Organism Comparison for EC {ec_number}', fontsize=16, fontweight='bold')

        # Extract data
        names = [c['organism'] for c in valid_data]
        avg_kms = [c.get('average_km', 0) for c in valid_data if c.get('average_km')]
        optimal_phs = [c.get('optimal_ph', 0) for c in valid_data if c.get('optimal_ph')]
        optimal_temps = [c.get('optimal_temperature', 0) for c in valid_data if c.get('optimal_temperature')]
        data_points = [c.get('data_points', 0) for c in valid_data]

        # Plot 1: Average Km comparison
        if avg_kms:
            ax1.bar(names, avg_kms)
            ax1.set_ylabel('Average Km (mM)')
            ax1.set_title('Average Km Comparison')
            ax1.tick_params(axis='x', rotation=45)

        # Plot 2: Optimal pH comparison
        if optimal_phs:
            ax2.bar(names, optimal_phs)
            ax2.set_ylabel('Optimal pH')
            ax2.set_title('Optimal pH Comparison')
            ax2.tick_params(axis='x', rotation=45)

        # Plot 3: Optimal temperature comparison
        if optimal_temps:
            ax3.bar(names, optimal_temps)
            ax3.set_ylabel('Optimal Temperature (°C)')
            ax3.set_title('Optimal Temperature Comparison')
            ax3.tick_params(axis='x', rotation=45)

        # Plot 4: Data points comparison
        ax4.bar(names, data_points)
        ax4.set_ylabel('Number of Data Points')
        ax4.set_title('Available Data Points')
        ax4.tick_params(axis='x', rotation=45)

        plt.tight_layout()

        # Save plot
        if save_path:
            plt.savefig(save_path, dpi=300, bbox_inches='tight')
            print(f"Organism comparison plot saved to {save_path}")

        if show_plot:
            plt.show()
        else:
            plt.close()

        return save_path or f"organism_comparison_{ec_number.replace('.', '_')}.png"

    except Exception as e:
        print(f"Error plotting organism comparison: {e}")
        return save_path


def plot_pH_profiles(ec_number: str, save_path: str = None, show_plot: bool = True) -> str:
    """Plot pH activity profiles for an enzyme."""
    validate_dependencies()

    try:
        # Get kinetic data
        km_data = get_km_values(ec_number)

        if not km_data:
            print(f"No pH data found for EC {ec_number}")
            return save_path

        # Parse data and extract pH information
        ph_kms = []
        ph_organisms = []

        for entry in km_data:
            parsed = parse_km_entry(entry)
            if 'ph' in parsed and 'km_value_numeric' in parsed:
                ph_kms.append((parsed['ph'], parsed['km_value_numeric']))
                ph_organisms.append(parsed.get('organism', 'Unknown'))

        if not ph_kms:
            print(f"No pH-Km data found for EC {ec_number}")
            return save_path

        # Create figure
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
        fig.suptitle(f'pH Activity Profiles for EC {ec_number}', fontsize=16, fontweight='bold')

        # Extract data
        ph_values = [item[0] for item in ph_kms]
        km_values = [item[1] for item in ph_kms]

        # Plot 1: pH vs Km scatter plot
        _scatter = ax1.scatter(ph_values, km_values, alpha=0.6, s=50)
        ax1.set_xlabel('pH')
        ax1.set_ylabel('Km (mM)')
        ax1.set_title('pH vs Km Values')
        ax1.grid(True, alpha=0.3)

        # Add trend line
        if len(ph_values) > 2:
            z = np.polyfit(ph_values, km_values, 1)
            p = np.poly1d(z)
            ax1.plot(ph_values, p(ph_values), "r--", alpha=0.8, label=f'Trend: y={z[0]:.3f}x+{z[1]:.3f}')
            ax1.legend()

        # Plot 2: pH distribution histogram
        ax2.hist(ph_values, bins=20, alpha=0.7, edgecolor='black')
        ax2.set_xlabel('pH')
        ax2.set_ylabel('Frequency')
        ax2.set_title('pH Distribution')
        ax2.axvline(np.mean(ph_values), color='red', linestyle='--', label=f'Mean: {np.mean(ph_values):.2f}')
        ax2.axvline(np.median(ph_values), color='blue', linestyle='--', label=f'Median: {np.median(ph_values):.2f}')
        ax2.legend()

        plt.tight_layout()

        # Save plot
        if save_path:
            plt.savefig(save_path, dpi=300, bbox_inches='tight')
            print(f"pH profile plot saved to {save_path}")

        if show_plot:
            plt.show()
        else:
            plt.close()

        return save_path or f"ph_profile_{ec_number.replace('.', '_')}.png"

    except Exception as e:
        print(f"Error plotting pH profiles: {e}")
        return save_path


def plot_temperature_profiles(ec_number: str, save_path: str = None, show_plot: bool = True) -> str:
    """Plot temperature activity profiles for an enzyme."""
    validate_dependencies()

    try:
        # Get kinetic data
        km_data = get_km_values(ec_number)

        if not km_data:
            print(f"No temperature data found for EC {ec_number}")
            return save_path

        # Parse data and extract temperature information
        temp_kms = []
        temp_organisms = []

        for entry in km_data:
            parsed = parse_km_entry(entry)
            if 'temperature' in parsed and 'km_value_numeric' in parsed:
                temp_kms.append((parsed['temperature'], parsed['km_value_numeric']))
                temp_organisms.append(parsed.get('organism', 'Unknown'))

        if not temp_kms:
            print(f"No temperature-Km data found for EC {ec_number}")
            return save_path

        # Create figure
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
        fig.suptitle(f'Temperature Activity Profiles for EC {ec_number}', fontsize=16, fontweight='bold')

        # Extract data
        temp_values = [item[0] for item in temp_kms]
        km_values = [item[1] for item in temp_kms]

        # Plot 1: Temperature vs Km scatter plot
        _scatter = ax1.scatter(temp_values, km_values, alpha=0.6, s=50)
        ax1.set_xlabel('Temperature (°C)')
        ax1.set_ylabel('Km (mM)')
        ax1.set_title('Temperature vs Km Values')
        ax1.grid(True, alpha=0.3)

        # Add trend line
        if len(temp_values) > 2:
            z = np.polyfit(temp_values, km_values, 2)  # Quadratic fit for temperature optima
            p = np.poly1d(z)
            x_smooth = np.linspace(min(temp_values), max(temp_values), 100)
            ax1.plot(x_smooth, p(x_smooth), "r--", alpha=0.8, label='Polynomial fit')

            # Find optimum temperature
            optimum_idx = np.argmin(p(x_smooth))
            optimum_temp = x_smooth[optimum_idx]
            ax1.axvline(optimum_temp, color='green', linestyle=':', label=f'Optimal: {optimum_temp:.1f}°C')
            ax1.legend()

        # Plot 2: Temperature distribution histogram
        ax2.hist(temp_values, bins=20, alpha=0.7, edgecolor='black')
        ax2.set_xlabel('Temperature (°C)')
        ax2.set_ylabel('Frequency')
        ax2.set_title('Temperature Distribution')
        ax2.axvline(np.mean(temp_values), color='red', linestyle='--', label=f'Mean: {np.mean(temp_values):.1f}°C')
        ax2.axvline(np.median(temp_values), color='blue', linestyle='--', label=f'Median: {np.median(temp_values):.1f}°C')
        ax2.legend()

        plt.tight_layout()

        # Save plot
        if save_path:
            plt.savefig(save_path, dpi=300, bbox_inches='tight')
            print(f"Temperature profile plot saved to {save_path}")

        if show_plot:
            plt.show()
        else:
            plt.close()

        return save_path or f"temperature_profile_{ec_number.replace('.', '_')}.png"

    except Exception as e:
        print(f"Error plotting temperature profiles: {e}")
        return save_path


def plot_substrate_specificity(ec_number: str, save_path: str = None, show_plot: bool = True) -> str:
    """Plot substrate specificity and affinity for an enzyme."""
    validate_dependencies()

    try:
        # Get substrate specificity data
        specificity = get_substrate_specificity(ec_number)

        if not specificity:
            print(f"No substrate specificity data found for EC {ec_number}")
            return save_path

        # Create figure
        fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
        fig.suptitle(f'Substrate Specificity for EC {ec_number}', fontsize=16, fontweight='bold')

        # Extract data
        substrates = [s['name'] for s in specificity]
        kms = [s['km'] for s in specificity if s.get('km')]
        data_points = [s['data_points'] for s in specificity]

        # Get top substrates for plotting
        if PANDAS_AVAILABLE and kms:
            df = pd.DataFrame({'Substrate': substrates, 'Km': kms, 'DataPoints': data_points})
            top_substrates = df.nlargest(15, 'DataPoints')  # Top 15 by data points

            # Plot 1: Km values for top substrates (sorted by affinity)
            top_sorted = top_substrates.sort_values('Km')
            ax1.barh(range(len(top_sorted)), top_sorted['Km'])
            ax1.set_yticks(range(len(top_sorted)))
            ax1.set_yticklabels([s[:30] + '...' if len(s) > 30 else s for s in top_sorted['Substrate']])
            ax1.set_xlabel('Km (mM)')
            ax1.set_title('Substrate Affinity (Lower Km = Higher Affinity)')
            ax1.invert_yaxis()  # Best affinity at top

            # Plot 2: Data points by substrate
            ax2.barh(range(len(top_sorted)), top_sorted['DataPoints'])
            ax2.set_yticks(range(len(top_sorted)))
            ax2.set_yticklabels([s[:30] + '...' if len(s) > 30 else s for s in top_sorted['Substrate']])
            ax2.set_xlabel('Number of Data Points')
            ax2.set_title('Data Availability by Substrate')
            ax2.invert_yaxis()

            # Plot 3: Km distribution
            ax3.hist(kms, bins=20, alpha=0.7, edgecolor='black')
            ax3.set_xlabel('Km (mM)')
            ax3.set_ylabel('Frequency')
            ax3.set_title('Km Value Distribution')
            ax3.axvline(np.mean(kms), color='red', linestyle='--', label=f'Mean: {np.mean(kms):.2f}')
            ax3.axvline(np.median(kms), color='blue', linestyle='--', label=f'Median: {np.median(kms):.2f}')
            ax3.legend()

            # Plot 4: Km vs Data Points scatter
            ax4.scatter(df['DataPoints'], df['Km'], alpha=0.6)
            ax4.set_xlabel('Number of Data Points')
            ax4.set_ylabel('Km (mM)')
            ax4.set_title('Km vs Data Points')
            ax4.grid(True, alpha=0.3)

        plt.tight_layout()

        # Save plot
        if save_path:
            plt.savefig(save_path, dpi=300, bbox_inches='tight')
            print(f"Substrate specificity plot saved to {save_path}")

        if show_plot:
            plt.show()
        else:
            plt.close()

        return save_path or f"substrate_specificity_{ec_number.replace('.', '_')}.png"

    except Exception as e:
        print(f"Error plotting substrate specificity: {e}")
        return save_path


def plot_michaelis_menten(ec_number: str, substrate: str = None, save_path: str = None, show_plot: bool = True) -> str:
    """Generate Michaelis-Menten curves for an enzyme."""
    validate_dependencies()

    try:
        # Get modeling parameters
        model_data = get_modeling_parameters(ec_number, substrate)

        if not model_data or model_data.get('error'):
            print(f"No modeling data found for EC {ec_number}")
            return save_path

        km = model_data.get('km')
        vmax = model_data.get('vmax')
        kcat = model_data.get('kcat')
        enzyme_conc = model_data.get('enzyme_conc', 1.0)

        if not km:
            print(f"No Km data available for plotting")
            return save_path

        # Create figure
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
        fig.suptitle(f'Michaelis-Menten Kinetics for EC {ec_number}' + (f' - {substrate}' if substrate else ''),
                     fontsize=16, fontweight='bold')

        # Generate substrate concentration range
        substrate_range = np.linspace(0, km * 5, 1000)

        # Calculate reaction rates
        if vmax:
            # Use actual Vmax if available
            rates = (vmax * substrate_range) / (km + substrate_range)
        elif kcat and enzyme_conc:
            # Calculate Vmax from kcat and enzyme concentration
            vmax_calc = kcat * enzyme_conc
            rates = (vmax_calc * substrate_range) / (km + substrate_range)
        else:
            # Use normalized Vmax = 1.0
            rates = substrate_range / (km + substrate_range)

        # Plot 1: Michaelis-Menten curve
        ax1.plot(substrate_range, rates, 'b-', linewidth=2, label='Michaelis-Menten')
        ax1.axhline(y=rates[-1] * 0.5, color='r', linestyle='--', alpha=0.7, label='0.5 × Vmax')
        ax1.axvline(x=km, color='g', linestyle='--', alpha=0.7, label=f'Km = {km:.2f}')
        ax1.set_xlabel('Substrate Concentration (mM)')
        ax1.set_ylabel('Reaction Rate')
        ax1.set_title('Michaelis-Menten Curve')
        ax1.legend()
        ax1.grid(True, alpha=0.3)

        # Add annotation for Km
        km_rate = (substrate_range[km == min(substrate_range, key=lambda x: abs(x-km))] *
                  (vmax if vmax else kcat * enzyme_conc if kcat else 1.0)) / (km +
                  substrate_range[km == min(substrate_range, key=lambda x: abs(x-km))])
        ax1.plot(km, km_rate, 'ro', markersize=8)

        # Plot 2: Lineweaver-Burk plot (double reciprocal)
        substrate_range_nonzero = substrate_range[substrate_range > 0]
        rates_nonzero = rates[substrate_range > 0]

        reciprocal_substrate = 1 / substrate_range_nonzero
        reciprocal_rate = 1 / rates_nonzero

        ax2.scatter(reciprocal_substrate, reciprocal_rate, alpha=0.6, s=10)

        # Fit linear regression
        z = np.polyfit(reciprocal_substrate, reciprocal_rate, 1)
        p = np.poly1d(z)
        x_fit = np.linspace(min(reciprocal_substrate), max(reciprocal_substrate), 100)
        ax2.plot(x_fit, p(x_fit), 'r-', linewidth=2, label=f'1/Vmax = {z[1]:.3f}')

        ax2.set_xlabel('1/[Substrate] (1/mM)')
        ax2.set_ylabel('1/Rate')
        ax2.set_title('Lineweaver-Burk Plot')
        ax2.legend()
        ax2.grid(True, alpha=0.3)

        # Add parameter information
        info_text = f"Km = {km:.3f} mM"
        if vmax:
            info_text += f"\nVmax = {vmax:.3f}"
        if kcat:
            info_text += f"\nkcat = {kcat:.3f} s⁻¹"
        if enzyme_conc:
            info_text += f"\n[Enzyme] = {enzyme_conc:.3f} μM"

        fig.text(0.02, 0.98, info_text, transform=fig.transFigure,
                fontsize=10, verticalalignment='top',
                bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))

        plt.tight_layout()

        # Save plot
        if save_path:
            plt.savefig(save_path, dpi=300, bbox_inches='tight')
            print(f"Michaelis-Menten plot saved to {save_path}")

        if show_plot:
            plt.show()
        else:
            plt.close()

        return save_path or f"michaelis_menten_{ec_number.replace('.', '_')}_{substrate or 'all'}.png"

    except Exception as e:
        print(f"Error plotting Michaelis-Menten: {e}")
        return save_path


def create_heatmap_data(ec_number: str, parameters: List[str] = None) -> Dict[str, Any]:
    """Create data for heatmap visualization."""
    validate_dependencies()

    try:
        # Get comparison data across organisms
        organisms = ["Escherichia coli", "Saccharomyces cerevisiae", "Bacillus subtilis",
                    "Homo sapiens", "Mus musculus", "Rattus norvegicus"]
        comparison = compare_across_organisms(ec_number, organisms)

        if not comparison:
            return None

        # Create heatmap data
        heatmap_data = {
            'organisms': [],
            'average_km': [],
            'optimal_ph': [],
            'optimal_temperature': [],
            'data_points': []
        }

        for comp in comparison:
            if comp.get('data_points', 0) > 0:
                heatmap_data['organisms'].append(comp['organism'])
                heatmap_data['average_km'].append(comp.get('average_km', 0))
                heatmap_data['optimal_ph'].append(comp.get('optimal_ph', 0))
                heatmap_data['optimal_temperature'].append(comp.get('optimal_temperature', 0))
                heatmap_data['data_points'].append(comp.get('data_points', 0))

        return heatmap_data

    except Exception as e:
        print(f"Error creating heatmap data: {e}")
        return None


def plot_heatmap(ec_number: str, save_path: str = None, show_plot: bool = True) -> str:
    """Create heatmap visualization of enzyme properties."""
    validate_dependencies()

    try:
        heatmap_data = create_heatmap_data(ec_number)

        if not heatmap_data or not heatmap_data['organisms']:
            print(f"No heatmap data available for EC {ec_number}")
            return save_path

        if not PANDAS_AVAILABLE:
            print("pandas required for heatmap plotting")
            return save_path

        # Create DataFrame for heatmap
        df = pd.DataFrame({
            'Organism': heatmap_data['organisms'],
            'Avg Km (mM)': heatmap_data['average_km'],
            'Optimal pH': heatmap_data['optimal_ph'],
            'Optimal Temp (°C)': heatmap_data['optimal_temperature'],
            'Data Points': heatmap_data['data_points']
        })

        # Normalize data for better visualization
        df_normalized = df.copy()
        for col in ['Avg Km (mM)', 'Optimal pH', 'Optimal Temp (°C)', 'Data Points']:
            if col in df.columns:
                df_normalized[col] = (df[col] - df[col].min()) / (df[col].max() - df[col].min())

        # Create figure
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))
        fig.suptitle(f'Enzyme Properties Heatmap for EC {ec_number}', fontsize=16, fontweight='bold')

        # Plot 1: Raw data heatmap
        heatmap_data_raw = df.set_index('Organism')[['Avg Km (mM)', 'Optimal pH', 'Optimal Temp (°C)', 'Data Points']].T
        sns.heatmap(heatmap_data_raw, annot=True, fmt='.2f', cmap='viridis', ax=ax1)
        ax1.set_title('Raw Values')

        # Plot 2: Normalized data heatmap
        heatmap_data_norm = df_normalized.set_index('Organism')[['Avg Km (mM)', 'Optimal pH', 'Optimal Temp (°C)', 'Data Points']].T
        sns.heatmap(heatmap_data_norm, annot=True, fmt='.2f', cmap='viridis', ax=ax2)
        ax2.set_title('Normalized Values (0-1)')

        plt.tight_layout()

        # Save plot
        if save_path:
            plt.savefig(save_path, dpi=300, bbox_inches='tight')
            print(f"Heatmap plot saved to {save_path}")

        if show_plot:
            plt.show()
        else:
            plt.close()

        return save_path or f"heatmap_{ec_number.replace('.', '_')}.png"

    except Exception as e:
        print(f"Error plotting heatmap: {e}")
        return save_path


def generate_summary_plots(ec_number: str, save_dir: str = None) -> List[str]:
    """Generate a comprehensive set of plots for an enzyme."""
    validate_dependencies()

    if save_dir is None:
        save_dir = f"enzyme_plots_{ec_number.replace('.', '_')}"

    # Create save directory
    Path(save_dir).mkdir(exist_ok=True)

    generated_files = []

    # Generate all plot types
    plot_functions = [
        ('kinetic_parameters', plot_kinetic_parameters),
        ('ph_profiles', plot_pH_profiles),
        ('temperature_profiles', plot_temperature_profiles),
        ('substrate_specificity', plot_substrate_specificity),
        ('heatmap', plot_heatmap),
    ]

    for plot_name, plot_func in plot_functions:
        try:
            save_path = f"{save_dir}/{plot_name}_{ec_number.replace('.', '_')}.png"
            result_path = plot_func(ec_number, save_path=save_path, show_plot=False)
            if result_path:
                generated_files.append(result_path)
                print(f"Generated {plot_name} plot")
            else:
                print(f"Failed to generate {plot_name} plot")
        except Exception as e:
            print(f"Error generating {plot_name} plot: {e}")

    # Generate organism comparison for common model organisms
    model_organisms = ["Escherichia coli", "Saccharomyces cerevisiae", "Homo sapiens"]
    try:
        save_path = f"{save_dir}/organism_comparison_{ec_number.replace('.', '_')}.png"
        result_path = plot_organism_comparison(ec_number, model_organisms, save_path=save_path, show_plot=False)
        if result_path:
            generated_files.append(result_path)
            print("Generated organism comparison plot")
    except Exception as e:
        print(f"Error generating organism comparison plot: {e}")

    # Generate Michaelis-Menten plot for most common substrate
    try:
        specificity = get_substrate_specificity(ec_number)
        if specificity:
            most_common = max(specificity, key=lambda x: x.get('data_points', 0))
            substrate_name = most_common['name'].split()[0]  # Take first word
            save_path = f"{save_dir}/michaelis_menten_{ec_number.replace('.', '_')}_{substrate_name}.png"
            result_path = plot_michaelis_menten(ec_number, substrate_name, save_path=save_path, show_plot=False)
            if result_path:
                generated_files.append(result_path)
                print(f"Generated Michaelis-Menten plot for {substrate_name}")
    except Exception as e:
        print(f"Error generating Michaelis-Menten plot: {e}")

    print(f"\nGenerated {len(generated_files)} plots in directory: {save_dir}")
    return generated_files


if __name__ == "__main__":
    # Example usage
    print("BRENDA Visualization Examples")
    print("=" * 40)

    try:
        ec_number = "1.1.1.1"  # Alcohol dehydrogenase

        print(f"\n1. Generating kinetic parameters plot for EC {ec_number}")
        plot_kinetic_parameters(ec_number, show_plot=False)

        print(f"\n2. Generating pH profile plot for EC {ec_number}")
        plot_pH_profiles(ec_number, show_plot=False)

        print(f"\n3. Generating substrate specificity plot for EC {ec_number}")
        plot_substrate_specificity(ec_number, show_plot=False)

        print(f"\n4. Generating Michaelis-Menten plot for EC {ec_number}")
        plot_michaelis_menten(ec_number, substrate="ethanol", show_plot=False)

        print(f"\n5. Generating organism comparison plot for EC {ec_number}")
        organisms = ["Escherichia coli", "Saccharomyces cerevisiae", "Homo sapiens"]
        plot_organism_comparison(ec_number, organisms, show_plot=False)

        print(f"\n6. Generating comprehensive summary plots for EC {ec_number}")
        summary_files = generate_summary_plots(ec_number, show_plot=False)
        print(f"Generated {len(summary_files)} summary plots")

    except Exception as e:
        print(f"Example failed: {e}")

"""
Enzyme Pathway Builder for Retrosynthetic Analysis

This module provides tools for constructing enzymatic pathways and
retrosynthetic trees using BRENDA database information.

Key features:
- Find enzymatic pathways for target products
- Build retrosynthetic trees from products
- Suggest enzyme substitutions and alternatives
- Calculate pathway feasibility and thermodynamics
- Optimize pathway conditions (pH, temperature, cofactors)
- Generate detailed pathway reports
- Support for metabolic engineering and synthetic biology

Installation:
    uv pip install networkx matplotlib pandas

Usage:
    from scripts.enzyme_pathway_builder import find_pathway_for_product, build_retrosynthetic_tree

    pathway = find_pathway_for_product("lactate", max_steps=3)
    tree = build_retrosynthetic_tree("lactate", depth=2)
"""

import re
import time
from typing import List, Dict, Any

try:
    import networkx as nx
    NETWORKX_AVAILABLE = True
except ImportError:
    print("Warning: networkx not installed. Install with: uv pip install networkx")
    NETWORKX_AVAILABLE = False

try:
    import pandas as pd  # noqa: F401  # availability probe for PANDAS_AVAILABLE
    PANDAS_AVAILABLE = True
except ImportError:
    print("Warning: pandas not installed. Install with: uv pip install pandas")
    PANDAS_AVAILABLE = False

try:
    import matplotlib.pyplot as plt
    MATPLOTLIB_AVAILABLE = True
except ImportError:
    print("Warning: matplotlib not installed. Install with: uv pip install matplotlib")
    MATPLOTLIB_AVAILABLE = False

try:
    from brenda_queries import (
        search_enzymes_by_product, search_enzymes_by_substrate,
        get_environmental_parameters, compare_across_organisms,
        search_by_pattern, get_cofactor_requirements,
        find_thermophilic_homologs, find_ph_stable_variants
    )
    BRENDA_QUERIES_AVAILABLE = True
except ImportError:
    print("Warning: brenda_queries not available")
    BRENDA_QUERIES_AVAILABLE = False


def validate_dependencies():
    """Validate that required dependencies are installed."""
    missing = []
    if not NETWORKX_AVAILABLE:
        missing.append("networkx")
    if not PANDAS_AVAILABLE:
        missing.append("pandas")
    if not BRENDA_QUERIES_AVAILABLE:
        missing.append("brenda_queries")
    if missing:
        raise ImportError(f"Missing required dependencies: {', '.join(missing)}")


# Common biochemical transformations with typical EC numbers
COMMON_TRANSFORMATIONS = {
    'oxidation': ['1.1.1'],      # Alcohol dehydrogenases
    'reduction': ['1.1.1'],      # Alcohol dehydrogenases
    'hydrolysis': ['3.1.1', '3.1.3'],  # Esterases, phosphatases
    'carboxylation': ['6.4.1'],   # Carboxylases
    'decarboxylation': ['4.1.1'], # Decarboxylases
    'transamination': ['2.6.1'],  # Aminotransferases
    'phosphorylation': ['2.7.1'], # Kinases
    'dephosphorylation': ['3.1.3'], # Phosphatases
    'isomerization': ['5.1.1', '5.3.1'], # Isomerases
    'ligation': ['6.3.1'],       # Ligases
    'transfer': ['2.1.1', '2.2.1', '2.4.1'], # Transferases
    'hydride_transfer': ['1.1.1', '1.2.1'],  # Oxidoreductases
    'group_transfer': ['2.1.1'],  # Methyltransferases
}

# Simple metabolite database (expanded for pathway building)
METABOLITE_DATABASE = {
    # Primary metabolites
    'glucose': {'formula': 'C6H12O6', 'mw': 180.16, 'class': 'sugar'},
    'fructose': {'formula': 'C6H12O6', 'mw': 180.16, 'class': 'sugar'},
    'galactose': {'formula': 'C6H12O6', 'mw': 180.16, 'class': 'sugar'},
    'pyruvate': {'formula': 'C3H4O3', 'mw': 90.08, 'class': 'carboxylic_acid'},
    'lactate': {'formula': 'C3H6O3', 'mw': 90.08, 'class': 'carboxylic_acid'},
    'acetate': {'formula': 'C2H4O2', 'mw': 60.05, 'class': 'carboxylic_acid'},
    'ethanol': {'formula': 'C2H6O', 'mw': 46.07, 'class': 'alcohol'},
    'acetaldehyde': {'formula': 'C2H4O', 'mw': 44.05, 'class': 'aldehyde'},
    'acetone': {'formula': 'C3H6O', 'mw': 58.08, 'class': 'ketone'},
    'glycerol': {'formula': 'C3H8O3', 'mw': 92.09, 'class': 'alcohol'},
    'ammonia': {'formula': 'NH3', 'mw': 17.03, 'class': 'inorganic'},
    'carbon dioxide': {'formula': 'CO2', 'mw': 44.01, 'class': 'inorganic'},
    'water': {'formula': 'H2O', 'mw': 18.02, 'class': 'inorganic'},
    'oxygen': {'formula': 'O2', 'mw': 32.00, 'class': 'inorganic'},
    'hydrogen': {'formula': 'H2', 'mw': 2.02, 'class': 'inorganic'},
    'nitrogen': {'formula': 'N2', 'mw': 28.01, 'class': 'inorganic'},
    'phosphate': {'formula': 'PO4', 'mw': 94.97, 'class': 'inorganic'},
    'sulfate': {'formula': 'SO4', 'mw': 96.06, 'class': 'inorganic'},

    # Amino acids
    'alanine': {'formula': 'C3H7NO2', 'mw': 89.09, 'class': 'amino_acid'},
    'glycine': {'formula': 'C2H5NO2', 'mw': 75.07, 'class': 'amino_acid'},
    'serine': {'formula': 'C3H7NO3', 'mw': 105.09, 'class': 'amino_acid'},
    'threonine': {'formula': 'C4H9NO3', 'mw': 119.12, 'class': 'amino_acid'},
    'aspartate': {'formula': 'C4H7NO4', 'mw': 133.10, 'class': 'amino_acid'},
    'glutamate': {'formula': 'C5H9NO4', 'mw': 147.13, 'class': 'amino_acid'},
    'asparagine': {'formula': 'C4H8N2O3', 'mw': 132.12, 'class': 'amino_acid'},
    'glutamine': {'formula': 'C5H10N2O3', 'mw': 146.15, 'class': 'amino_acid'},
    'lysine': {'formula': 'C6H14N2O2', 'mw': 146.19, 'class': 'amino_acid'},
    'arginine': {'formula': 'C6H14N4O2', 'mw': 174.20, 'class': 'amino_acid'},
    'histidine': {'formula': 'C6H9N3O2', 'mw': 155.16, 'class': 'amino_acid'},
    'phenylalanine': {'formula': 'C9H11NO2', 'mw': 165.19, 'class': 'amino_acid'},
    'tyrosine': {'formula': 'C9H11NO3', 'mw': 181.19, 'class': 'amino_acid'},
    'tryptophan': {'formula': 'C11H12N2O2', 'mw': 204.23, 'class': 'amino_acid'},
    'leucine': {'formula': 'C6H13NO2', 'mw': 131.18, 'class': 'amino_acid'},
    'isoleucine': {'formula': 'C6H13NO2', 'mw': 131.18, 'class': 'amino_acid'},
    'valine': {'formula': 'C5H11NO2', 'mw': 117.15, 'class': 'amino_acid'},
    'methionine': {'formula': 'C5H11NO2S', 'mw': 149.21, 'class': 'amino_acid'},
    'cysteine': {'formula': 'C3H7NO2S', 'mw': 121.16, 'class': 'amino_acid'},
    'proline': {'formula': 'C5H9NO2', 'mw': 115.13, 'class': 'amino_acid'},

    # Nucleotides (simplified)
    'atp': {'formula': 'C10H16N5O13P3', 'mw': 507.18, 'class': 'nucleotide'},
    'adp': {'formula': 'C10H15N5O10P2', 'mw': 427.20, 'class': 'nucleotide'},
    'amp': {'formula': 'C10H14N5O7P', 'mw': 347.22, 'class': 'nucleotide'},
    'nad': {'formula': 'C21H27N7O14P2', 'mw': 663.43, 'class': 'cofactor'},
    'nadh': {'formula': 'C21H29N7O14P2', 'mw': 665.44, 'class': 'cofactor'},
    'nadp': {'formula': 'C21H28N7O17P3', 'mw': 743.44, 'class': 'cofactor'},
    'nadph': {'formula': 'C21H30N7O17P3', 'mw': 745.45, 'class': 'cofactor'},
    'fadh2': {'formula': 'C21H30N7O14P2', 'mw': 785.55, 'class': 'cofactor'},
    'fadx': {'formula': 'C21H20N4O2', 'mw': 350.36, 'class': 'cofactor'},

    # Common organic acids
    'malate': {'formula': 'C4H6O5', 'mw': 134.09, 'class': 'carboxylic_acid'},
    'oxaloacetate': {'formula': 'C4H4O5', 'mw': 132.07, 'class': 'carboxylic_acid'},
    'succinate': {'formula': 'C4H6O4', 'mw': 118.09, 'class': 'carboxylic_acid'},
    'fumarate': {'formula': 'C4H4O4', 'mw': 116.07, 'class': 'carboxylic_acid'},
    'oxalosuccinate': {'formula': 'C6H6O7', 'mw': 190.12, 'class': 'carboxylic_acid'},
    'alpha-ketoglutarate': {'formula': 'C5H6O5', 'mw': 146.11, 'class': 'carboxylic_acid'},

    # Energy carriers
    'acetyl-coa': {'formula': 'C23H38N7O17P3S', 'mw': 809.51, 'class': 'cofactor'},
    'coenzyme-a': {'formula': 'C21H36N7O16P3S', 'mw': 767.54, 'class': 'cofactor'},
}

# Common cofactors and their roles
COFACTOR_ROLES = {
    'nad+': {'role': 'oxidation', 'oxidation_state': '+1'},
    'nadh': {'role': 'reduction', 'oxidation_state': '0'},
    'nadp+': {'role': 'oxidation', 'oxidation_state': '+1'},
    'nadph': {'role': 'reduction', 'oxidation_state': '0'},
    'fadx': {'role': 'oxidation', 'oxidation_state': '0'},
    'fadh2': {'role': 'reduction', 'oxidation_state': '-2'},
    'atp': {'role': 'phosphorylation', 'oxidation_state': '0'},
    'adp': {'role': 'energy', 'oxidation_state': '0'},
    'amp': {'role': 'energy', 'oxidation_state': '0'},
    'acetyl-coa': {'role': 'acetylation', 'oxidation_state': '0'},
    'coenzyme-a': {'role': 'thiolation', 'oxidation_state': '0'},
}


def identify_metabolite(metabolite_name: str) -> Dict[str, Any]:
    """Identify a metabolite from the database or create entry."""
    metabolite_name = metabolite_name.lower().strip()

    # Check if it's in the database
    if metabolite_name in METABOLITE_DATABASE:
        return {'name': metabolite_name, **METABOLITE_DATABASE[metabolite_name]}

    # Simple formula extraction from common patterns
    formula_patterns = {
        r'c(\d+)h(\d+)o(\d+)': lambda m: f"C{m[0]}H{m[1]}O{m[2]}",
        r'c(\d+)h(\d+)n(\d+)o(\d+)': lambda m: f"C{m[0]}H{m[1]}N{m[2]}O{m[3]}",
    }

    for pattern, formatter in formula_patterns.items():
        match = re.search(pattern, metabolite_name)
        if match:
            formula = formatter(match.groups())
            # Estimate molecular weight (C=12, H=1, N=14, O=16)
            mw = 0
            elements = re.findall(r'([A-Z])(\d*)', formula)
            for elem, count in elements:
                count = int(count) if count else 1
                if elem == 'C':
                    mw += count * 12.01
                elif elem == 'H':
                    mw += count * 1.008
                elif elem == 'N':
                    mw += count * 14.01
                elif elem == 'O':
                    mw += count * 16.00
                elif elem == 'P':
                    mw += count * 30.97
                elif elem == 'S':
                    mw += count * 32.07

            return {
                'name': metabolite_name,
                'formula': formula,
                'mw': mw,
                'class': 'unknown'
            }

    # Fallback - unknown metabolite
    return {
        'name': metabolite_name,
        'formula': 'Unknown',
        'mw': 0,
        'class': 'unknown'
    }


def infer_transformation_type(substrate: str, product: str) -> List[str]:
    """Infer the type of transformation based on substrate and product."""
    substrate_info = identify_metabolite(substrate)
    product_info = identify_metabolite(product)

    transformations = []

    # Check for oxidation/reduction patterns
    if 'alcohol' in substrate_info.get('class', '') and 'carboxylic_acid' in product_info.get('class', ''):
        transformations.append('oxidation')
    elif 'aldehyde' in substrate_info.get('class', '') and 'alcohol' in product_info.get('class', ''):
        transformations.append('reduction')
    elif 'alcohol' in substrate_info.get('class', '') and 'aldehyde' in product_info.get('class', ''):
        transformations.append('oxidation')

    # Check for phosphorylation/dephosphorylation
    if 'phosphate' in product and 'phosphate' not in substrate:
        transformations.append('phosphorylation')
    elif 'phosphate' in substrate and 'phosphate' not in product:
        transformations.append('dephosphorylation')

    # Check for carboxylation/decarboxylation
    if 'co2' in product and 'co2' not in substrate:
        transformations.append('carboxylation')
    elif 'co2' in substrate and 'co2' not in product:
        transformations.append('decarboxylation')

    # Check for hydrolysis (simple heuristic)
    if 'ester' in substrate.lower() and ('carboxylic_acid' in product_info.get('class', '') or 'alcohol' in product_info.get('class', '')):
        transformations.append('hydrolysis')

    # Check for transamination
    if 'amino_acid' in product_info.get('class', '') and 'amino_acid' not in substrate_info.get('class', ''):
        transformations.append('transamination')

    # Default to generic transformation
    if not transformations:
        transformations.append('generic')

    return transformations


def find_enzymes_for_transformation(substrate: str, product: str, limit: int = 10) -> List[Dict[str, Any]]:
    """Find enzymes that catalyze a specific transformation."""
    validate_dependencies()

    # Infer transformation types
    transformations = infer_transformation_type(substrate, product)

    all_enzymes = []

    # Try to find enzymes by product
    try:
        product_enzymes = search_enzymes_by_product(product, limit=limit)
        for enzyme in product_enzymes:
            # Check if substrate is in the reactants
            if substrate.lower() in enzyme.get('reaction', '').lower():
                enzyme['transformation'] = transformations[0] if transformations else 'generic'
                enzyme['substrate'] = substrate
                enzyme['product'] = product
                enzyme['confidence'] = 'high'
                all_enzymes.append(enzyme)
        time.sleep(0.5)  # Rate limiting
    except Exception as e:
        print(f"Error searching enzymes by product: {e}")

    # Try to find enzymes by substrate
    try:
        substrate_enzymes = search_enzymes_by_substrate(substrate, limit=limit)
        for enzyme in substrate_enzymes:
            # Check if product is mentioned in substrate data (limited approach)
            enzyme['transformation'] = transformations[0] if transformations else 'generic'
            enzyme['substrate'] = substrate
            enzyme['product'] = product
            enzyme['confidence'] = 'medium'
            all_enzymes.append(enzyme)
        time.sleep(0.5)  # Rate limiting
    except Exception as e:
        print(f"Error searching enzymes by substrate: {e}")

    # If no enzymes found, try common EC numbers for transformation types
    if not all_enzymes and transformations:
        for trans_type in transformations:
            if trans_type in COMMON_TRANSFORMATIONS:
                for ec_prefix in COMMON_TRANSFORMATIONS[trans_type]:
                    # This is a simplified approach - in practice you'd want
                    # to query the specific EC numbers with more detail
                    try:
                        generic_enzymes = search_by_pattern(trans_type, limit=5)
                        for enzyme in generic_enzymes:
                            enzyme['transformation'] = trans_type
                            enzyme['substrate'] = substrate
                            enzyme['product'] = product
                            enzyme['confidence'] = 'low'
                            all_enzymes.append(enzyme)
                        time.sleep(0.5)
                        break
                    except Exception as e:
                        print(f"Error searching for transformation type {trans_type}: {e}")

    # Remove duplicates and sort by confidence
    unique_enzymes = []
    seen = set()
    for enzyme in all_enzymes:
        key = (enzyme.get('ec_number', ''), enzyme.get('organism', ''))
        if key not in seen:
            seen.add(key)
            unique_enzymes.append(enzyme)

    # Sort by confidence (high > medium > low)
    confidence_order = {'high': 3, 'medium': 2, 'low': 1}
    unique_enzymes.sort(key=lambda x: confidence_order.get(x.get('confidence', 'low'), 0), reverse=True)

    return unique_enzymes[:limit]


def find_pathway_for_product(product: str, max_steps: int = 3, starting_materials: List[str] = None) -> Dict[str, Any]:
    """Find enzymatic pathways to synthesize a target product."""
    validate_dependencies()

    if starting_materials is None:
        # Common starting materials
        starting_materials = ['glucose', 'pyruvate', 'acetate', 'ethanol', 'glycerol']

    pathway = {
        'target': product,
        'max_steps': max_steps,
        'starting_materials': starting_materials,
        'steps': [],
        'alternative_pathways': [],
        'warnings': [],
        'confidence': 0
    }

    # Simple breadth-first search for pathway
    from collections import deque

    queue = deque([(product, 0, [product])])  # (current_metabolite, step_count, pathway)
    visited = set()

    while queue and len(pathway['steps']) == 0:
        current_metabolite, step_count, current_path = queue.popleft()

        if current_metabolite in visited or step_count >= max_steps:
            continue

        visited.add(current_metabolite)

        # Check if current metabolite is a starting material
        if current_metabolite.lower() in [sm.lower() for sm in starting_materials]:
            # Found a complete pathway
            pathway['steps'] = []
            for i in range(len(current_path) - 1):
                substrate = current_path[i + 1]
                product_step = current_path[i]
                enzymes = find_enzymes_for_transformation(substrate, product_step, limit=5)

                if enzymes:
                    pathway['steps'].append({
                        'step_number': i + 1,
                        'substrate': substrate,
                        'product': product_step,
                        'enzymes': enzymes,
                        'transformation': infer_transformation_type(substrate, product_step)
                    })
                else:
                    pathway['warnings'].append(f"No enzymes found for step: {substrate} -> {product_step}")

            pathway['confidence'] = 0.8  # High confidence for found pathway
            break

        # Try to find enzymes that produce current metabolite
        if step_count < max_steps:
            # Generate possible substrates (simplified - in practice you'd need metabolic knowledge)
            possible_substrates = []

            # Try common metabolic precursors
            common_precursors = ['glucose', 'pyruvate', 'acetate', 'ethanol', 'acetyl-CoA', 'oxaloacetate']
            for precursor in common_precursors:
                enzymes = find_enzymes_for_transformation(precursor, current_metabolite, limit=2)
                if enzymes:
                    possible_substrates.append(precursor)
                    pathway['alternative_pathways'].append({
                        'precursor': precursor,
                        'product': current_metabolite,
                        'enzymes': enzymes
                    })

            # Add found substrates to queue
            for substrate in possible_substrates:
                if substrate not in current_path:
                    new_path = [substrate] + current_path
                    queue.append((substrate, step_count + 1, new_path))

        time.sleep(0.2)  # Rate limiting

    # If no complete pathway found, create partial pathway
    if not pathway['steps'] and pathway['alternative_pathways']:
        # Create best guess pathway from alternatives
        best_alternative = max(pathway['alternative_pathways'],
                               key=lambda x: len(x.get('enzymes', [])))
        pathway['steps'] = [{
            'step_number': 1,
            'substrate': best_alternative['precursor'],
            'product': best_alternative['product'],
            'enzymes': best_alternative['enzymes'],
            'transformation': infer_transformation_type(best_alternative['precursor'], best_alternative['product'])
        }]
        pathway['confidence'] = 0.3  # Low confidence for partial pathway
        pathway['warnings'].append("Partial pathway only - complete synthesis route not found")

    elif not pathway['steps']:
        pathway['warnings'].append("No enzymatic pathway found for target product")
        pathway['confidence'] = 0.1

    return pathway


def build_retrosynthetic_tree(target: str, depth: int = 2) -> Dict[str, Any]:
    """Build a retrosynthetic tree for a target molecule."""
    validate_dependencies()

    tree = {
        'target': target,
        'depth': depth,
        'nodes': {target: {'level': 0, 'children': [], 'enzymes': []}},
        'edges': [],
        'alternative_routes': []
    }

    # Build tree recursively
    def build_node_recursive(metabolite: str, current_depth: int, parent: str = None) -> None:
        if current_depth >= depth:
            return

        # Find enzymes that can produce this metabolite
        potential_precursors = ['glucose', 'pyruvate', 'acetate', 'ethanol', 'acetyl-CoA',
                                'oxaloacetate', 'alpha-ketoglutarate', 'malate']

        for precursor in potential_precursors:
            enzymes = find_enzymes_for_transformation(precursor, metabolite, limit=3)

            if enzymes:
                # Add precursor as node if not exists
                if precursor not in tree['nodes']:
                    tree['nodes'][precursor] = {
                        'level': current_depth + 1,
                        'children': [],
                        'enzymes': enzymes
                    }
                    tree['nodes'][metabolite]['children'].append(precursor)
                    tree['edges'].append({
                        'from': precursor,
                        'to': metabolite,
                        'enzymes': enzymes,
                        'transformation': infer_transformation_type(precursor, metabolite)
                    })

                # Recursively build tree
                if current_depth + 1 < depth:
                    build_node_recursive(precursor, current_depth + 1, metabolite)

        # Try common metabolic transformations
        if current_depth < depth - 1:
            transformations = ['oxidation', 'reduction', 'hydrolysis', 'carboxylation', 'decarboxylation']
            for trans in transformations:
                try:
                    generic_enzymes = search_by_pattern(trans, limit=2)
                    if generic_enzymes:
                        # Create hypothetical precursor
                        hypothetical_precursor = f"precursor_{trans}_{metabolite}"
                        tree['nodes'][hypothetical_precursor] = {
                            'level': current_depth + 1,
                            'children': [],
                            'enzymes': generic_enzymes,
                            'hypothetical': True
                        }
                        tree['nodes'][metabolite]['children'].append(hypothetical_precursor)
                        tree['edges'].append({
                            'from': hypothetical_precursor,
                            'to': metabolite,
                            'enzymes': generic_enzymes,
                            'transformation': trans,
                            'hypothetical': True
                        })
                except Exception as e:
                    print(f"Error in retrosynthetic search for {trans}: {e}")

        time.sleep(0.3)  # Rate limiting

    # Start building from target
    build_node_recursive(target, 0)

    # Calculate tree statistics
    tree['total_nodes'] = len(tree['nodes'])
    tree['total_edges'] = len(tree['edges'])
    tree['max_depth'] = max(node['level'] for node in tree['nodes'].values()) if tree['nodes'] else 0

    return tree


def suggest_enzyme_substitutions(ec_number: str, criteria: Dict[str, Any] = None) -> List[Dict[str, Any]]:
    """Suggest alternative enzymes with improved properties."""
    validate_dependencies()

    if criteria is None:
        criteria = {
            'min_temperature': 30,
            'max_temperature': 70,
            'min_ph': 6.0,
            'max_ph': 8.0,
            'min_thermostability': 40,
            'prefer_organisms': ['Escherichia coli', 'Saccharomyces cerevisiae', 'Bacillus subtilis']
        }

    substitutions = []

    # Get organisms for the target enzyme
    try:
        organisms = compare_across_organisms(ec_number, criteria['prefer_organisms'])
        time.sleep(0.5)
    except Exception as e:
        print(f"Error comparing organisms: {e}")
        organisms = []

    # Find thermophilic homologs if temperature is a criterion
    if criteria.get('min_thermostability'):
        try:
            thermophilic = find_thermophilic_homologs(ec_number, criteria['min_thermostability'])
            time.sleep(0.5)

            for enzyme in thermophilic:
                enzyme['substitution_reason'] = f"Thermostable (optimal temp: {enzyme['optimal_temperature']}°C)"
                enzyme['score'] = 8.0 if enzyme['optimal_temperature'] >= criteria['min_thermostability'] else 6.0
                substitutions.append(enzyme)
        except Exception as e:
            print(f"Error finding thermophilic homologs: {e}")

    # Find pH-stable variants
    if criteria.get('min_ph') or criteria.get('max_ph'):
        try:
            ph_stable = find_ph_stable_variants(ec_number, criteria.get('min_ph'), criteria.get('max_ph'))
            time.sleep(0.5)

            for enzyme in ph_stable:
                enzyme['substitution_reason'] = f"pH stable ({enzyme['stability_type']} range: {enzyme['ph_range']})"
                enzyme['score'] = 7.5
                substitutions.append(enzyme)
        except Exception as e:
            print(f"Error finding pH-stable variants: {e}")

    # Add organism comparison results
    for org_data in organisms:
        if org_data.get('data_points', 0) > 0:
            org_data['substitution_reason'] = f"Well-characterized in {org_data['organism']}"
            org_data['score'] = 6.5 if org_data['organism'] in criteria['prefer_organisms'] else 5.0
            substitutions.append(org_data)

    # Sort by score
    substitutions.sort(key=lambda x: x.get('score', 0), reverse=True)

    return substitutions[:10]  # Return top 10 suggestions


def calculate_pathway_feasibility(pathway: Dict[str, Any]) -> Dict[str, Any]:
    """Calculate feasibility scores and potential issues for a pathway."""
    validate_dependencies()

    feasibility = {
        'overall_score': 0,
        'step_scores': [],
        'warnings': [],
        'recommendations': [],
        'thermodynamic_feasibility': 0,
        'enzyme_availability': 0,
        'cofactor_requirements': [],
        'optimal_conditions': {}
    }

    if not pathway.get('steps'):
        feasibility['warnings'].append("No steps in pathway")
        feasibility['overall_score'] = 0.1
        return feasibility

    total_score = 0
    step_scores = []

    for step in pathway['steps']:
        step_score = 0
        enzymes = step.get('enzymes', [])

        # Score based on number of available enzymes
        if len(enzymes) >= 3:
            step_score += 3  # Multiple enzyme options
        elif len(enzymes) >= 1:
            step_score += 2  # At least one enzyme
        else:
            step_score += 0  # No enzymes
            feasibility['warnings'].append(f"No enzymes found for step: {step['substrate']} -> {step['product']}")

        # Score based on enzyme confidence
        if enzymes:
            high_confidence = sum(1 for e in enzymes if e.get('confidence') == 'high')
            confidence_bonus = min(high_confidence, 2)  # Max 2 points for confidence
            step_score += confidence_bonus

        # Check for industrial viability
        industrial_organisms = ['Escherichia coli', 'Saccharomyces cerevisiae', 'Bacillus subtilis']
        industrial_enzymes = sum(1 for e in enzymes if e.get('organism') in industrial_organisms)
        if industrial_enzymes > 0:
            step_score += 1

        # Cap step score at 5
        step_score = min(step_score, 5)
        step_scores.append(step_score)
        total_score += step_score

        # Analyze cofactor requirements
        try:
            for enzyme in enzymes:
                ec_number = enzyme.get('ec_number', '')
                if ec_number:
                    cofactors = get_cofactor_requirements(ec_number)
                    for cofactor in cofactors:
                        if cofactor['name'] not in [c['name'] for c in feasibility['cofactor_requirements']]:
                            feasibility['cofactor_requirements'].append(cofactor)
            time.sleep(0.3)
        except Exception as e:
            print(f"Error analyzing cofactors: {e}")

    feasibility['step_scores'] = step_scores
    feasibility['enzyme_availability'] = total_score / (len(step_scores) * 5)  # Normalize to 0-1
    feasibility['overall_score'] = feasibility['enzyme_availability'] * 0.7  # Weight enzyme availability

    # Thermodynamic feasibility (simplified heuristic)
    pathway_length = len(pathway['steps'])
    if pathway_length <= 2:
        feasibility['thermodynamic_feasibility'] = 0.8  # Short pathways are often feasible
    elif pathway_length <= 4:
        feasibility['thermodynamic_feasibility'] = 0.6
    else:
        feasibility['thermodynamic_feasibility'] = 0.4  # Long pathways may have thermodynamic issues

    # Overall feasibility is weighted combination
    feasibility['overall_score'] = (
        feasibility['enzyme_availability'] * 0.6 +
        feasibility['thermodynamic_feasibility'] * 0.4
    )

    # Generate recommendations
    if feasibility['overall_score'] < 0.3:
        feasibility['warnings'].append("Low overall pathway feasibility")
        feasibility['recommendations'].append("Consider alternative starting materials or target molecules")
    elif feasibility['overall_score'] < 0.6:
        feasibility['warnings'].append("Moderate pathway feasibility")
        feasibility['recommendations'].append("Consider enzyme engineering or cofactor recycling")

    if feasibility['cofactor_requirements']:
        feasibility['recommendations'].append("Implement cofactor recycling system for: " +
                                            ", ".join([c['name'] for c in feasibility['cofactor_requirements']]))

    return feasibility


def optimize_pathway_conditions(pathway: Dict[str, Any]) -> Dict[str, Any]:
    """Suggest optimal conditions for the entire pathway."""
    validate_dependencies()

    optimization = {
        'optimal_temperature': 30.0,  # Default
        'optimal_ph': 7.0,           # Default
        'temperature_range': (20, 40),  # Default
        'ph_range': (6.5, 7.5),         # Default
        'cofactor_system': [],
        'organism_compatibility': {},
        'process_recommendations': []
    }

    temperatures = []
    phs = []
    organism_preferences = {}

    # Collect environmental data from all enzymes
    for step in pathway.get('steps', []):
        for enzyme in step.get('enzymes', []):
            ec_number = enzyme.get('ec_number', '')
            organism = enzyme.get('organism', '')

            if ec_number:
                try:
                    env_params = get_environmental_parameters(ec_number)
                    time.sleep(0.3)

                    if env_params.get('optimal_temperature'):
                        temperatures.append(env_params['optimal_temperature'])
                    if env_params.get('optimal_ph'):
                        phs.append(env_params['optimal_ph'])

                    # Track organism preferences
                    if organism not in organism_preferences:
                        organism_preferences[organism] = {
                            'temperature_optima': [],
                            'ph_optima': [],
                            'step_count': 0
                        }

                    organism_preferences[organism]['step_count'] += 1
                    if env_params.get('optimal_temperature'):
                        organism_preferences[organism]['temperature_optima'].append(env_params['optimal_temperature'])
                    if env_params.get('optimal_ph'):
                        organism_preferences[organism]['ph_optima'].append(env_params['optimal_ph'])

                except Exception as e:
                    print(f"Error getting environmental parameters for {ec_number}: {e}")

    # Calculate optimal conditions
    if temperatures:
        optimization['optimal_temperature'] = sum(temperatures) / len(temperatures)
        optimization['temperature_range'] = (min(temperatures) - 5, max(temperatures) + 5)

    if phs:
        optimization['optimal_ph'] = sum(phs) / len(phs)
        optimization['ph_range'] = (min(phs) - 0.5, max(phs) + 0.5)

    # Find best organism compatibility
    for organism, data in organism_preferences.items():
        if data['temperature_optima'] and data['ph_optima']:
            organism_preferences[organism]['avg_temp'] = sum(data['temperature_optima']) / len(data['temperature_optima'])
            organism_preferences[organism]['avg_ph'] = sum(data['ph_optima']) / len(data['ph_optima'])
            organism_preferences[organism]['compatibility_score'] = data['step_count']

    # Sort organisms by compatibility
    compatible_organisms = sorted(
        [(org, data) for org, data in organism_preferences.items() if data.get('compatibility_score', 0) > 0],
        key=lambda x: x[1]['compatibility_score'],
        reverse=True
    )

    optimization['organism_compatibility'] = dict(compatible_organisms[:5])  # Top 5 organisms

    # Generate process recommendations
    if len(optimization['organism_compatibility']) > 1:
        optimization['process_recommendations'].append("Consider multi-organism system or enzyme cocktails")

    if optimization['temperature_range'][1] - optimization['temperature_range'][0] > 30:
        optimization['process_recommendations'].append("Consider temperature gradient or staged process")

    if optimization['ph_range'][1] - optimization['ph_range'][0] > 2:
        optimization['process_recommendations'].append("Consider pH control system or buffer optimization")

    # Cofactor system optimization
    cofactor_types = {}
    for step in pathway.get('steps', []):
        for enzyme in step.get('enzymes', []):
            ec_number = enzyme.get('ec_number', '')
            if ec_number:
                try:
                    cofactors = get_cofactor_requirements(ec_number)
                    for cofactor in cofactors:
                        cofactor_type = cofactor.get('type', 'other')
                        if cofactor_type not in cofactor_types:
                            cofactor_types[cofactor_type] = []
                        if cofactor['name'] not in cofactor_types[cofactor_type]:
                            cofactor_types[cofactor_type].append(cofactor['name'])
                    time.sleep(0.3)
                except Exception as e:
                    print(f"Error getting cofactors for {ec_number}: {e}")

    optimization['cofactor_system'] = cofactor_types

    return optimization


def generate_pathway_report(pathway: Dict[str, Any], filename: str = None) -> str:
    """Generate a comprehensive pathway report."""
    validate_dependencies()

    if filename is None:
        target_name = pathway.get('target', 'pathway').replace(' ', '_').lower()
        filename = f"pathway_report_{target_name}.txt"

    # Calculate feasibility and optimization
    feasibility = calculate_pathway_feasibility(pathway)
    optimization = optimize_pathway_conditions(pathway)

    report = []
    report.append("=" * 80)
    report.append(f"ENZYMATIC PATHWAY REPORT")
    report.append("=" * 80)

    # Overview
    report.append(f"\nTARGET PRODUCT: {pathway.get('target', 'Unknown')}")
    report.append(f"PATHWAY LENGTH: {len(pathway.get('steps', []))} steps")
    report.append(f"OVERALL FEASIBILITY: {feasibility['overall_score']:.2f}/1.00")

    # Pathway steps
    if pathway.get('steps'):
        report.append("\n" + "=" * 40)
        report.append("PATHWAY STEPS")
        report.append("=" * 40)

        for i, step in enumerate(pathway['steps'], 1):
            report.append(f"\nStep {i}: {step['substrate']} -> {step['product']}")
            report.append(f"Transformation: {', '.join(step.get('transformation', ['Unknown']))}")

            if step.get('enzymes'):
                report.append(f"Available enzymes: {len(step['enzymes'])}")
                for j, enzyme in enumerate(step['enzymes'][:3], 1):  # Top 3 enzymes
                    report.append(f"  {j}. EC {enzyme.get('ec_number', 'Unknown')} - {enzyme.get('organism', 'Unknown')}")
                    report.append(f"     Confidence: {enzyme.get('confidence', 'Unknown')}")
                    if enzyme.get('reaction'):
                        report.append(f"     Reaction: {enzyme['reaction'][:100]}...")

                if len(step['enzymes']) > 3:
                    report.append(f"  ... and {len(step['enzymes']) - 3} additional enzymes")
            else:
                report.append("  No enzymes found for this step")

            if feasibility.get('step_scores') and i-1 < len(feasibility['step_scores']):
                report.append(f"Step feasibility score: {feasibility['step_scores'][i-1]}/5.0")

    # Cofactor requirements
    if feasibility.get('cofactor_requirements'):
        report.append("\n" + "=" * 40)
        report.append("COFACTOR REQUIREMENTS")
        report.append("=" * 40)

        for cofactor in feasibility['cofactor_requirements']:
            report.append(f"- {cofactor['name']} ({cofactor.get('type', 'Unknown')})")
            report.append(f"  Organism: {cofactor.get('organism', 'Unknown')}")
            report.append(f"  EC Number: {cofactor.get('ec_number', 'Unknown')}")

    # Optimal conditions
    report.append("\n" + "=" * 40)
    report.append("OPTIMAL CONDITIONS")
    report.append("=" * 40)

    report.append(f"Temperature: {optimization['optimal_temperature']:.1f}°C")
    report.append(f"pH: {optimization['optimal_ph']:.1f}")
    report.append(f"Temperature range: {optimization['temperature_range'][0]:.1f} - {optimization['temperature_range'][1]:.1f}°C")
    report.append(f"pH range: {optimization['ph_range'][0]:.1f} - {optimization['ph_range'][1]:.1f}")

    if optimization.get('organism_compatibility'):
        report.append("\nCompatible organisms (by preference):")
        for organism, data in list(optimization['organism_compatibility'].items())[:3]:
            report.append(f"- {organism} (compatibility score: {data.get('compatibility_score', 0)})")
            if data.get('avg_temp'):
                report.append(f"  Optimal temperature: {data['avg_temp']:.1f}°C")
            if data.get('avg_ph'):
                report.append(f"  Optimal pH: {data['avg_ph']:.1f}")

    # Warnings and recommendations
    if feasibility.get('warnings'):
        report.append("\n" + "=" * 40)
        report.append("WARNINGS")
        report.append("=" * 40)

        for warning in feasibility['warnings']:
            report.append(f"⚠️  {warning}")

    if feasibility.get('recommendations'):
        report.append("\n" + "=" * 40)
        report.append("RECOMMENDATIONS")
        report.append("=" * 40)

        for rec in feasibility['recommendations']:
            report.append(f"💡 {rec}")

    if optimization.get('process_recommendations'):
        for rec in optimization['process_recommendations']:
            report.append(f"🔧 {rec}")

    # Alternative pathways
    if pathway.get('alternative_pathways'):
        report.append("\n" + "=" * 40)
        report.append("ALTERNATIVE ROUTES")
        report.append("=" * 40)

        for alt in pathway['alternative_pathways'][:5]:  # Top 5 alternatives
            report.append(f"\n{alt['precursor']} -> {alt['product']}")
            report.append(f"Enzymes available: {len(alt.get('enzymes', []))}")
            for enzyme in alt.get('enzymes', [])[:2]:  # Top 2 enzymes
                report.append(f"  - {enzyme.get('ec_number', 'Unknown')} ({enzyme.get('organism', 'Unknown')})")

    # Feasibility analysis
    report.append("\n" + "=" * 40)
    report.append("FEASIBILITY ANALYSIS")
    report.append("=" * 40)

    report.append(f"Enzyme availability score: {feasibility['enzyme_availability']:.2f}/1.00")
    report.append(f"Thermodynamic feasibility: {feasibility['thermodynamic_feasibility']:.2f}/1.00")

    # Write report to file
    with open(filename, 'w') as f:
        f.write('\n'.join(report))

    print(f"Pathway report saved to {filename}")
    return filename


def visualize_pathway(pathway: Dict[str, Any], save_path: str = None) -> str:
    """Create a visual representation of the pathway."""
    validate_dependencies()

    if not NETWORKX_AVAILABLE or not MATPLOTLIB_AVAILABLE:
        print("networkx and matplotlib required for pathway visualization")
        return save_path or "pathway_visualization.png"

    try:
        # Create directed graph
        G = nx.DiGraph()

        # Add nodes and edges
        for step in pathway.get('steps', []):
            substrate = step['substrate']
            product = step['product']
            enzymes = step.get('enzymes', [])

            G.add_node(substrate, type='substrate')
            G.add_node(product, type='product')

            # Add edge with enzyme information
            edge_label = f"{len(enzymes)} enzymes"
            if enzymes:
                primary_ec = enzymes[0].get('ec_number', 'Unknown')
                edge_label += f"\nEC {primary_ec}"

            G.add_edge(substrate, product, label=edge_label)

        # Create figure
        plt.figure(figsize=(12, 8))

        # Layout
        pos = nx.spring_layout(G, k=2, iterations=50)

        # Draw nodes
        substrate_nodes = [n for n, d in G.nodes(data=True) if d.get('type') == 'substrate']
        product_nodes = [n for n, d in G.nodes(data=True) if d.get('type') == 'product']
        intermediate_nodes = [n for n in G.nodes() if n not in substrate_nodes and n not in product_nodes]

        nx.draw_networkx_nodes(G, pos, nodelist=substrate_nodes, node_color='lightblue', node_size=1500)
        nx.draw_networkx_nodes(G, pos, nodelist=product_nodes, node_color='lightgreen', node_size=1500)
        nx.draw_networkx_nodes(G, pos, nodelist=intermediate_nodes, node_color='lightyellow', node_size=1200)

        # Draw edges
        nx.draw_networkx_edges(G, pos, edge_color='gray', arrows=True, arrowsize=20)

        # Draw labels
        nx.draw_networkx_labels(G, pos, font_size=10, font_weight='bold')

        # Draw edge labels
        edge_labels = nx.get_edge_attributes(G, 'label')
        nx.draw_networkx_edge_labels(G, pos, edge_labels, font_size=8)

        # Add title
        plt.title(f"Enzymatic Pathway to {pathway.get('target', 'Target')}", fontsize=14, fontweight='bold')

        # Add legend
        plt.scatter([], [], c='lightblue', s=150, label='Starting Materials')
        plt.scatter([], [], c='lightyellow', s=120, label='Intermediates')
        plt.scatter([], [], c='lightgreen', s=150, label='Products')
        plt.legend()

        plt.axis('off')
        plt.tight_layout()

        # Save or show
        if save_path:
            plt.savefig(save_path, dpi=300, bbox_inches='tight')
            print(f"Pathway visualization saved to {save_path}")
        else:
            plt.show()

        plt.close()
        return save_path or "pathway_visualization.png"

    except Exception as e:
        print(f"Error visualizing pathway: {e}")
        return save_path or "pathway_visualization.png"


if __name__ == "__main__":
    # Example usage
    print("Enzyme Pathway Builder Examples")
    print("=" * 50)

    try:
        # Example 1: Find pathway for lactate
        print("\n1. Finding pathway for lactate production:")
        pathway = find_pathway_for_product("lactate", max_steps=3)
        print(f"Found pathway with {len(pathway['steps'])} steps")
        print(f"Feasibility: {pathway['confidence']:.2f}")

        # Example 2: Build retrosynthetic tree
        print("\n2. Building retrosynthetic tree for ethanol:")
        tree = build_retrosynthetic_tree("ethanol", depth=2)
        print(f"Tree has {tree['total_nodes']} nodes and {tree['total_edges']} edges")

        # Example 3: Suggest enzyme substitutions
        print("\n3. Suggesting enzyme substitutions for alcohol dehydrogenase:")
        substitutions = suggest_enzyme_substitutions("1.1.1.1")
        for sub in substitutions[:3]:
            print(f"  - {sub.get('organism', 'Unknown')}: {sub.get('substitution_reason', 'No reason')}")

        # Example 4: Calculate feasibility
        print("\n4. Calculating pathway feasibility:")
        feasibility = calculate_pathway_feasibility(pathway)
        print(f"Overall score: {feasibility['overall_score']:.2f}")
        print(f"Warnings: {len(feasibility['warnings'])}")

        # Example 5: Generate pathway report
        print("\n5. Generating pathway report:")
        report_file = generate_pathway_report(pathway)
        print(f"Report saved to: {report_file}")

        # Example 6: Visualize pathway
        print("\n6. Visualizing pathway:")
        viz_file = visualize_pathway(pathway, "example_pathway.png")
        print(f"Visualization saved to: {viz_file}")

    except Exception as e:
        print(f"Example failed: {e}")