{
  "skill": "alterlab-bioservices",
  "evals": [
    {
      "id": "uniprot-to-kegg-mapping",
      "prompt": "I have the UniProt accession P43403 (ZAP70). I want to map it to its KEGG gene ID and then pull every KEGG pathway that gene participates in. Can you script this?",
      "expected_output": "Invokes alterlab-bioservices. Uses the UniProt service mapping() (fr='UniProtKB_AC-ID', to='KEGG', query='P43403') to convert the accession, then uses the KEGG service (get_pathway_by_gene / lookfor_pathway) to enumerate pathways for the gene. This is a cross-database workflow spanning UniProt and KEGG, which is bioservices' core strength.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "mapping" },
        { "type": "behavior", "value": "Chains UniProt identifier mapping into KEGG pathway lookup rather than treating it as a single-database call." }
      ]
    },
    {
      "id": "compound-cross-reference",
      "prompt": "Find the KEGG compound ID for Geldanamycin and cross-reference it to its ChEBI identifier so I can pull the chemical structure later.",
      "expected_output": "Invokes alterlab-bioservices. Uses KEGG.find('compound', 'Geldanamycin') to get cpd:C11222, then KEGG.conv('chebi', 'compound') to map KEGG -> ChEBI. Notes that there is no direct bioservices KEGG -> ChEMBL convenience method and that a ChEMBL ID would require a separate route. Multi-database compound cross-referencing is in scope.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "ChEBI" }
      ]
    },
    {
      "id": "ppi-string-query",
      "prompt": "I want the reported protein-protein interaction partners of ZAP70 in human. How do I query that programmatically?",
      "expected_output": "Invokes alterlab-bioservices. Uses the STRING service — STRING().get_interaction_partners('ZAP70', species=9606), or get_interactions for edges within a set — and notes that the old PSICQUIC and BioGRID wrappers were removed in bioservices 1.14, so PSICQUIC().query(...) no longer imports.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "STRING" }
      ]
    },
    {
      "id": "go-term-annotation",
      "prompt": "Pull the GO term definition for GO:0003824 and the GO annotations for UniProt P43403 in one go using QuickGO.",
      "expected_output": "Invokes alterlab-bioservices. Uses the QuickGO service: g.get_go_terms('GO:0003824') for the term and g.Annotation(geneProductId='UniProtKB:P43403', includeFields='goName', limit=100) for annotations, reading the 'results' list of the returned dict. Gene Ontology queries via QuickGO are an in-scope bioservices service.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "QuickGO" }
      ]
    },
    {
      "id": "near-miss-gget",
      "prompt": "I just need a quick one-off lookup of the human BRCA1 gene record from Ensembl — nothing fancy, no cross-database mapping.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-gget. The ask is a single quick single-database lookup with no multi-service integration or identifier mapping, which the lightweight gget tool handles. bioservices is reserved for workflows that hit multiple databases together or map identifiers across services.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-gget" }
      ]
    },
    {
      "id": "near-miss-biopython",
      "prompt": "I already downloaded a FASTA file of protein sequences. Parse it, translate the ORFs, and write out a multi-record GenBank file.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-biopython. The ask is local sequence parsing, translation, and file-format conversion with no web-service or database query, which is biopython's domain. bioservices is for hitting remote bioinformatics web services.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-biopython" }
      ]
    }
  ]
}

# BioServices: Identifier Mapping Guide

This document provides comprehensive information about converting identifiers between different biological databases using BioServices.

## Table of Contents

1. [Overview](#overview)
2. [UniProt Mapping Service](#uniprot-mapping-service)
3. [Compound Identifier Mapping](#compound-identifier-mapping)
4. [KEGG Identifier Conversions](#kegg-identifier-conversions)
5. [Common Mapping Patterns](#common-mapping-patterns)
6. [Troubleshooting](#troubleshooting)

---

## Overview

Biological databases use different identifier systems. Cross-referencing requires mapping between these systems. BioServices provides multiple approaches:

1. **UniProt Mapping**: Comprehensive protein/gene ID conversion (`UniProt.mapping`)
2. **KEGG**: `KEGG.conv` conversions plus cross-references embedded in entries
3. **Compounds**: KEGG -> ChEBI via `KEGG.conv`; KEGG/ChEBI -> ChEMBL needs an
   external route (see [Compound Identifier Mapping](#compound-identifier-mapping))

Note: the old `bioservices.UniChem` compound-mapping helpers were dropped in 2022,
and the EBI PICR service was retired years ago — neither is usable through
bioservices today.

---

## UniProt Mapping Service

The UniProt mapping service is the most comprehensive tool for protein and gene identifier conversion.

### Basic Usage

```python
from bioservices import UniProt

u = UniProt()

# Map single ID
job = u.mapping(
    fr="UniProtKB_AC-ID",    # Source database
    to="KEGG",                # Target database
    query="P43403"            # Identifier to convert
)

print(job)
# {'results': [{'from': 'P43403', 'to': 'hsa:7535'}], 'failedIds': []}
```

> **Return shape (bioservices >= 1.10).** `mapping()` hands back UniProt's own
> ID-mapping job payload: a dict with `results` (a list of `{"from": ..., "to": ...}`
> rows) and `failedIds`. It is **not** keyed by source ID. Older code that does
> `result[uniprot_id]` raises `KeyError`. Collapse it once and reuse:
>
> ```python
> from collections import defaultdict
>
> def as_dict(job):
>     """{'P43403': ['hsa:7535'], ...} from a UniProt mapping payload."""
>     out = defaultdict(list)
>     for row in (job or {}).get("results", []):
>         out[row["from"]].append(row["to"])
>     return dict(out)
> ```
>
> For some targets (e.g. `to="UniProtKB"`) each `to` value is a full entry object rather
> than a plain accession — check the type before treating it as a string.

### Batch Mapping

```python
# Map multiple IDs (a list or a comma-separated string both work)
ids = ["P43403", "P04637", "P53779"]
job = u.mapping(fr="UniProtKB_AC-ID", to="KEGG", query=ids)

for uniprot_id, kegg_ids in as_dict(job).items():   # as_dict from the note above
    print(f"{uniprot_id} → {kegg_ids}")

print("unmapped:", job["failedIds"])
```

### Supported Database Pairs

UniProt supports mapping between 100+ database pairs. Key ones include:

#### Protein/Gene Databases

| Source Format | Code | Target Format | Code |
|---------------|------|---------------|------|
| UniProtKB AC/ID | `UniProtKB_AC-ID` | KEGG | `KEGG` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | Ensembl | `Ensembl` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | Ensembl Protein | `Ensembl_Protein` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | Ensembl Transcript | `Ensembl_Transcript` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | RefSeq Protein | `RefSeq_Protein` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | RefSeq Nucleotide | `RefSeq_Nucleotide` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | GeneID (Entrez) | `GeneID` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | HGNC | `HGNC` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | MGI | `MGI` |
| KEGG | `KEGG` | UniProtKB | `UniProtKB` |
| Ensembl | `Ensembl` | UniProtKB | `UniProtKB` |
| GeneID | `GeneID` | UniProtKB | `UniProtKB` |

#### Structural Databases

| Source | Code | Target | Code |
|--------|------|--------|------|
| UniProtKB AC/ID | `UniProtKB_AC-ID` | PDB | `PDB` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | Pfam | `Pfam` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | InterPro | `InterPro` |
| PDB | `PDB` | UniProtKB | `UniProtKB` |

#### Expression & Proteomics

| Source | Code | Target | Code |
|--------|------|--------|------|
| UniProtKB AC/ID | `UniProtKB_AC-ID` | PRIDE | `PRIDE` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | ProteomicsDB | `ProteomicsDB` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | PaxDb | `PaxDb` |

#### Organism-Specific

| Source | Code | Target | Code |
|--------|------|--------|------|
| UniProtKB AC/ID | `UniProtKB_AC-ID` | FlyBase | `FlyBase` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | WormBase | `WormBase` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | SGD | `SGD` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | ZFIN | `ZFIN` |

#### Other Useful Mappings

| Source | Code | Target | Code |
|--------|------|--------|------|
| UniProtKB AC/ID | `UniProtKB_AC-ID` | GO | `GO` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | Reactome | `Reactome` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | STRING | `STRING` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | BioGRID | `BioGRID` |
| UniProtKB AC/ID | `UniProtKB_AC-ID` | OMA | `OMA` |

### Complete List of Database Codes

To get the complete, up-to-date list:

```python
from bioservices import UniProt

u = UniProt()

# This information is in the UniProt REST API documentation
# Common patterns:
# - Source databases typically end in source database name
# - UniProtKB uses "UniProtKB_AC-ID" or "UniProtKB"
# - Most other databases use their standard abbreviation
```

### Common Database Codes Reference

**Gene/Protein Identifiers:**
- `UniProtKB_AC-ID`: UniProt accession/ID
- `UniProtKB`: UniProt accession
- `KEGG`: KEGG gene IDs (e.g., hsa:7535)
- `GeneID`: NCBI Gene (Entrez) IDs
- `Ensembl`: Ensembl gene IDs
- `Ensembl_Protein`: Ensembl protein IDs
- `Ensembl_Transcript`: Ensembl transcript IDs
- `RefSeq_Protein`: RefSeq protein IDs (NP_)
- `RefSeq_Nucleotide`: RefSeq nucleotide IDs (NM_)

**Gene Nomenclature:**
- `HGNC`: Human Gene Nomenclature Committee
- `MGI`: Mouse Genome Informatics
- `RGD`: Rat Genome Database
- `SGD`: Saccharomyces Genome Database
- `FlyBase`: Drosophila database
- `WormBase`: C. elegans database
- `ZFIN`: Zebrafish database

**Structure:**
- `PDB`: Protein Data Bank
- `Pfam`: Protein families
- `InterPro`: Protein domains
- `SUPFAM`: Superfamily
- `PROSITE`: Protein motifs

**Pathways & Networks:**
- `Reactome`: Reactome pathways
- `BioCyc`: BioCyc pathways
- `PathwayCommons`: Pathway Commons
- `STRING`: Protein-protein networks
- `BioGRID`: Interaction database

### Mapping Examples

#### UniProt → KEGG

```python
from bioservices import UniProt

u = UniProt()

# Single mapping
job = u.mapping(fr="UniProtKB_AC-ID", to="KEGG", query="P43403")
print(job)  # {'results': [{'from': 'P43403', 'to': 'hsa:7535'}], 'failedIds': []}
```

#### KEGG → UniProt

```python
# Reverse mapping. Non-UniProt sources map only *to* UniProtKB /
# UniProtKB-Swiss-Prot / UniParc, so "to='UniProtKB_AC-ID'" is rejected here.
job = u.mapping(fr="KEGG", to="UniProtKB", query="hsa:7535")
print([r["from"] for r in job["results"]])  # ['hsa:7535']
```

#### UniProt → Ensembl

```python
# To Ensembl gene IDs
job = u.mapping(fr="UniProtKB_AC-ID", to="Ensembl", query="P43403")
print(as_dict(job))  # {'P43403': ['ENSG00000115085.x']}

# To Ensembl protein IDs
job = u.mapping(fr="UniProtKB_AC-ID", to="Ensembl_Protein", query="P43403")
print(as_dict(job))  # {'P43403': ['ENSP00000381359.x']}
```

#### UniProt → PDB

```python
# Find 3D structures
job = u.mapping(fr="UniProtKB_AC-ID", to="PDB", query="P04637")
print(as_dict(job))  # {'P04637': ['1A1U', '1AIE', '1C26', ...]}
```

#### UniProt → RefSeq

```python
# Get RefSeq protein IDs
job = u.mapping(fr="UniProtKB_AC-ID", to="RefSeq_Protein", query="P43403")
print(as_dict(job))  # {'P43403': ['NP_001070.2']}
```

#### Gene Name → UniProt (via search, then mapping)

```python
# First search for gene
search_result = u.search("gene:ZAP70 AND organism_id:9606", frmt="tsv", columns="accession")
lines = search_result.strip().split("\n")
if len(lines) > 1:
    uniprot_id = lines[1].split("\t")[0]

    # Then map to other databases (as_dict collapses the job payload)
    kegg_ids = as_dict(u.mapping(fr="UniProtKB_AC-ID", to="KEGG", query=uniprot_id))
    print(kegg_ids)
```

---

## Compound Identifier Mapping

> **Heads-up:** bioservices used to expose UniChem compound-mapping helpers
> (`UniChem.get_compound_id_from_kegg`, `get_src_compound_ids`,
> `get_all_compound_ids`). **These were dropped from bioservices in 2022 and no
> longer exist.** There is also no KEGG -> ChEMBL route in the KEGG Web Service.
> The patterns below use what actually works today.

### KEGG -> ChEBI (supported via bioservices)

```python
from bioservices import KEGG

k = KEGG()

# conv returns a dict keyed by the source IDs, e.g. {'cpd:C11222': 'chebi:5292', ...}
conv = k.conv("chebi", "compound")
chebi_id = conv.get("cpd:C11222")  # -> 'chebi:5292' (Geldanamycin)

# ChEBI IDs are also embedded in the KEGG compound entry text:
entry = k.get("cpd:C11222")
for line in entry.split("\n"):
    if "ChEBI:" in line:
        chebi_id = line.split("ChEBI:")[1].strip().split()[0]  # '5292'
        break
```

### KEGG / ChEBI -> ChEMBL (needs an external route)

bioservices cannot do this mapping. Obtain the ChEMBL ID with a separate tool:

- **ChEMBL web service / `chembl_webresource_client`** (the official ChEMBL client), or
- the **live UniChem REST API** at `https://www.ebi.ac.uk/unichem/` queried directly
  (e.g. via `requests`).

UniChem source-database IDs (for direct REST calls; verify against current UniChem
docs as the set evolves): ChEMBL=1, DrugBank=2, PDB=3, KEGG=6, ChEBI=7, PubChem=22.

The live UniChem REST API ( `https://www.ebi.ac.uk/unichem/` ) returns all source
IDs for a given compound from one request; consult its current OpenAPI docs for the
exact endpoint and parameter shape before wiring it up, since the schema has changed
over time. The `chembl_webresource_client` package (installed separately from
bioservices) is the other supported route into ChEMBL.

---

## KEGG Identifier Conversions

KEGG entries contain cross-references that can be extracted by parsing.

### Extract Database Links from KEGG Entry

```python
from bioservices import KEGG

k = KEGG()

# Get compound entry
entry = k.get("cpd:C11222")

# Parse for specific database
chebi_id = None
uniprot_ids = []

for line in entry.split("\n"):
    if "ChEBI:" in line:
        # Extract ChEBI ID
        parts = line.split("ChEBI:")
        if len(parts) > 1:
            chebi_id = parts[1].strip().split()[0]

# For genes/proteins
gene_entry = k.get("hsa:7535")
for line in gene_entry.split("\n"):
    if line.startswith("            "):  # Database links section
        if "UniProt:" in line:
            parts = line.split("UniProt:")
            if len(parts) > 1:
                uniprot_id = parts[1].strip()
                uniprot_ids.append(uniprot_id)
```

### KEGG Gene ID Components

KEGG gene IDs have format `organism:gene_id`:

```python
kegg_id = "hsa:7535"
organism, gene_id = kegg_id.split(":")

print(f"Organism: {organism}")  # hsa (human)
print(f"Gene ID: {gene_id}")    # 7535
```

### KEGG Pathway to Genes

```python
k = KEGG()

# Get pathway entry
pathway = k.get("path:hsa04660")

# Parse for gene list
genes = []
in_gene_section = False

for line in pathway.split("\n"):
    if line.startswith("GENE"):
        in_gene_section = True

    if in_gene_section:
        if line.startswith(" " * 12):  # Gene line
            parts = line.strip().split()
            if parts:
                gene_id = parts[0]
                genes.append(f"hsa:{gene_id}")
        elif not line.startswith(" "):
            break

print(f"Found {len(genes)} genes")
```

---

## Common Mapping Patterns

### Pattern 1: Gene Symbol → Multiple Database IDs

```python
from bioservices import UniProt

def gene_symbol_to_ids(gene_symbol, organism="9606"):
    """Convert gene symbol to multiple database IDs."""
    u = UniProt()

    # Search for gene
    query = f"gene:{gene_symbol} AND organism_id:{organism}"
    result = u.search(query, frmt="tsv", columns="accession")

    lines = result.strip().split("\n")
    if len(lines) < 2:
        return None

    uniprot_id = lines[1].split("\t")[0]

    # Map to multiple databases
    def targets(to_db):
        job = u.mapping(fr="UniProtKB_AC-ID", to=to_db, query=uniprot_id)
        return [row["to"] for row in job.get("results", [])]

    ids = {
        'uniprot': uniprot_id,
        'kegg': targets("KEGG"),
        'ensembl': targets("Ensembl"),
        'refseq': targets("RefSeq_Protein"),
        'pdb': targets("PDB"),
    }

    return ids

# Usage
ids = gene_symbol_to_ids("ZAP70")
print(ids)
```

### Pattern 2: Compound Name → KEGG + ChEBI IDs

```python
from bioservices import KEGG

def compound_name_to_ids(compound_name):
    """Search a compound by name and return its KEGG and ChEBI IDs.

    Note: KEGG -> ChEMBL has no bioservices route; obtain a ChEMBL ID
    separately via chembl_webresource_client or the live UniChem REST API.
    """
    k = KEGG()

    # Search KEGG
    results = k.find("compound", compound_name)
    if not results:
        return None

    # Extract KEGG ID
    kegg_id = results.strip().split("\n")[0].split("\t")[0].replace("cpd:", "")

    # KEGG -> ChEBI via conv (supported), with entry-parse fallback
    chebi_id = k.conv("chebi", "compound").get(f"cpd:{kegg_id}")
    if not chebi_id:
        entry = k.get(f"cpd:{kegg_id}")
        for line in entry.split("\n"):
            if "ChEBI:" in line:
                chebi_id = line.split("ChEBI:")[1].strip().split()[0]
                break

    return {
        'kegg': kegg_id,
        'chebi': chebi_id,
        'chembl': None,  # not resolvable via bioservices; see note above
    }

# Usage
ids = compound_name_to_ids("Geldanamycin")
print(ids)
```

### Pattern 3: Batch ID Conversion with Error Handling

```python
from bioservices import UniProt

def safe_batch_mapping(ids, from_db, to_db, chunk_size=100):
    """Map IDs in chunks, returning {source_id: [target_ids]}."""
    u = UniProt()
    all_results = {}

    for i in range(0, len(ids), chunk_size):
        chunk = ids[i:i + chunk_size]

        try:
            job = u.mapping(fr=from_db, to=to_db, query=",".join(chunk))
            for src, targets in as_dict(job).items():   # as_dict defined above
                all_results.setdefault(src, []).extend(targets)
            print(f"✓ Processed {min(i + chunk_size, len(ids))}/{len(ids)}")

        except Exception as e:
            print(f"✗ Error at chunk {i}: {e}")

            # Retry the chunk one ID at a time so one bad ID doesn't lose the rest
            for single_id in chunk:
                try:
                    job = u.mapping(fr=from_db, to=to_db, query=single_id)
                    all_results.update(as_dict(job))
                except Exception as inner:
                    print(f"  ✗ {single_id}: {inner}")
                    all_results.setdefault(single_id, [])

    return all_results

# Usage
uniprot_ids = ["P43403", "P04637", "P53779", "INVALID123"]
mapping = safe_batch_mapping(uniprot_ids, "UniProtKB_AC-ID", "KEGG")
```

### Pattern 4: Multi-Hop Mapping

Sometimes you need to map through intermediate databases:

```python
from bioservices import UniProt

def multi_hop_mapping(gene_symbol, organism="9606"):
    """Gene symbol → UniProt → KEGG → Pathways."""
    u = UniProt()
    k = KEGG()

    # Step 1: Gene symbol → UniProt
    query = f"gene:{gene_symbol} AND organism_id:{organism}"
    result = u.search(query, frmt="tsv", columns="accession")

    lines = result.strip().split("\n")
    if len(lines) < 2:
        return None

    uniprot_id = lines[1].split("\t")[0]

    # Step 2: UniProt → KEGG
    job = u.mapping(fr="UniProtKB_AC-ID", to="KEGG", query=uniprot_id)
    kegg_ids = [row["to"] for row in job.get("results", []) if row["from"] == uniprot_id]
    if not kegg_ids:
        return None

    kegg_id = kegg_ids[0]

    # Step 3: KEGG → Pathways
    organism_code, gene_id = kegg_id.split(":")
    pathways = k.get_pathway_by_gene(gene_id, organism_code)

    return {
        'gene': gene_symbol,
        'uniprot': uniprot_id,
        'kegg': kegg_id,
        'pathways': pathways
    }

# Usage
result = multi_hop_mapping("TP53")
print(result)
```

---

## Troubleshooting

### Issue 1: No Mapping Found

**Symptom:** Mapping returns empty or None

**Solutions:**
1. Verify source ID exists in source database
2. Check database code spelling
3. Try reverse mapping
4. Some IDs may not have mappings in all databases

```python
job = u.mapping(fr="UniProtKB_AC-ID", to="KEGG", query="P43403")

if not job or not job.get("results"):
    print("No mapping found. Try:")
    print("1. Verify the ID exists: u.search('accession:P43403')")
    print("2. Check whether the protein has a KEGG annotation")
    print("3. Inspect job['failedIds'] for IDs UniProt rejected outright")
```

### Issue 2: Too Many IDs in Batch

**Symptom:** Batch mapping fails or times out

**Solution:** Split into smaller chunks

```python
def chunked_mapping(ids, from_db, to_db, chunk_size=50):
    all_results = {}

    for i in range(0, len(ids), chunk_size):
        chunk = ids[i:i + chunk_size]
        job = u.mapping(fr=from_db, to=to_db, query=",".join(chunk))
        all_results.update(as_dict(job))

    return all_results
```

### Issue 3: Multiple Target IDs

**Symptom:** One source ID maps to multiple target IDs

**Solution:** Handle as list

```python
job = u.mapping(fr="UniProtKB_AC-ID", to="PDB", query="P04637")
# as_dict(job) -> {'P04637': ['1A1U', '1AIE', '1C26', ...]}

pdb_ids = as_dict(job)['P04637']
print(f"Found {len(pdb_ids)} PDB structures")

for pdb_id in pdb_ids:
    print(f"  {pdb_id}")
```

### Issue 4: Organism Ambiguity

**Symptom:** Gene symbol maps to multiple organisms

**Solution:** Always specify organism in searches

```python
# Bad: Ambiguous
result = u.search("gene:TP53")  # Many organisms have TP53

# Good: Specific
result = u.search("gene:TP53 AND organism:9606")  # Human only
```

### Issue 5: Deprecated IDs

**Symptom:** Old database IDs don't map

**Solution:** Update to current IDs first

```python
# Check if ID is current
entry = u.retrieve("P43403", frmt="txt")

# Look for secondary accessions
for line in entry.split("\n"):
    if line.startswith("AC"):
        print(line)  # Shows primary and secondary accessions
```

---

## Best Practices

1. **Always validate inputs** before batch processing
2. **Handle None/empty results** gracefully
3. **Use chunking** for large ID lists (50-100 per chunk)
4. **Cache results** for repeated queries
5. **Specify organism** when possible to avoid ambiguity
6. **Log failures** in batch processing for later retry
7. **Add delays** between large batches to respect API limits

```python
import time

def polite_batch_mapping(ids, from_db, to_db):
    """Batch mapping with rate limiting."""
    results = {}

    for i in range(0, len(ids), 50):
        chunk = ids[i:i + 50]
        job = u.mapping(fr=from_db, to=to_db, query=",".join(chunk))
        results.update(as_dict(job))

        time.sleep(0.5)  # Be nice to the API

    return results
```

---

For complete working examples, see:
- `scripts/batch_id_converter.py`: Command-line batch conversion tool
- `workflow_patterns.md`: Integration into larger workflows

# BioServices: Complete Services Reference

This document provides a comprehensive reference for all major services available in BioServices, including key methods, parameters, and use cases.

## Protein & Gene Resources

### UniProt

Protein sequence and functional information database.

**Initialization:**
```python
from bioservices import UniProt
u = UniProt(verbose=False)
```

**Key Methods:**

- `search(query, frmt="tsv", columns=None, include_isoforms=False, sort="score", limit=None, size=25, database="uniprotkb", progress=False)`
  - Search UniProt with flexible query syntax
  - `frmt`: `"tsv"` (default), `"xlsx"`, `"fasta"`, `"json"`, `"gff"`. The pre-2022
    `"tab"`, `"txt"`, `"xml"` and `"rdf"` values are no longer accepted here.
  - `columns`: comma-separated **UniProt return-field names**, e.g.
    `"accession,gene_names,organism_name,length"` (the old display labels such as
    `"id,genes,organism"` were replaced in the June-2022 API)
  - Returns: string in the requested format (all pages concatenated)

- `retrieve(uniprot_id, frmt="json", database="uniprot", include=False)`
  - Retrieve specific UniProt entry (or a list of them)
  - `frmt`: "json" (default), "txt", "xml", "rdf", "gff", "fasta"
  - Returns: entry data in the requested format; a list when given a list of IDs

- `mapping(fr="UniProtKB_AC-ID", to="KEGG", query="P43403", polling_interval_seconds=3, max_waiting_time=100)`
  - Submit an ID-mapping job, poll it, and return the assembled payload
  - `query`: single ID, comma-separated string, or list
  - Returns: `{"results": [{"from": src, "to": tgt}, ...], "failedIds": [...]}` —
    **not** a `{source: [targets]}` dict; build that yourself if you need it
  - `valid_mapping` (property) lists the legal `fr`/`to` pairs straight from UniProt

- `get_df(entries, nChunk=100, organism=None, limit=10, columns=None)`
  - Build a pandas DataFrame for a list of accessions

**Common `columns` values:** accession, id, gene_names, gene_primary, organism_name,
organism_id, protein_name, length, sequence, mass, go_id, ec, xref_pdb, xref_kegg

**Use cases:**
- Protein sequence retrieval for BLAST
- Functional annotation lookup
- Cross-database identifier mapping
- Batch protein information retrieval

---

### KEGG (Kyoto Encyclopedia of Genes and Genomes)

Metabolic pathways, genes, and organisms database.

**Initialization:**
```python
from bioservices import KEGG
k = KEGG()
k.organism = "hsa"  # Set default organism
```

**Key Methods:**

- `list(database)`
  - List entries in KEGG database
  - `database`: "organism", "pathway", "module", "disease", "drug", "compound"
  - Returns: Multi-line string with entries

- `find(database, query)`
  - Search database by keywords
  - Returns: List of matching entries with IDs

- `get(entry_id)`
  - Retrieve entry by ID
  - Supports genes, pathways, compounds, etc.
  - Returns: Raw entry text

- `parse(data)`
  - Parse KEGG entry into dictionary
  - Returns: Dict with structured data

- `lookfor_organism(name)`
  - Search organisms by name pattern
  - Returns: List of matching organism codes

- `lookfor_pathway(name)`
  - Search pathways by name
  - Returns: List of pathway IDs

- `get_pathway_by_gene(gene_id, organism)`
  - Find pathways containing a gene, e.g. `get_pathway_by_gene("7535", "hsa")`
  - Returns the parsed PATHWAY block — a **dict** `{pathway_id: pathway_name}`
    (`{"hsa04064": "NF-kappa B signaling pathway", ...}`), despite the docstring's
    "list of pathway Ids". Iterating it yields the IDs; `.items()` gives you the
    names for free. IDs come back without the `path:` prefix.

- `parse_kgml_pathway(pathway_id)`
  - Parse pathway KGML for interactions
  - Returns: Dict with "entries" and "relations"

- `pathway2sif(pathway_id)`
  - Extract Simple Interaction Format data
  - Filters for activation/inhibition
  - Returns: List of interaction tuples

**Organism codes:**
- hsa: Homo sapiens
- mmu: Mus musculus
- dme: Drosophila melanogaster
- sce: Saccharomyces cerevisiae
- eco: Escherichia coli

**Use cases:**
- Pathway analysis and visualization
- Gene function annotation
- Metabolic network reconstruction
- Protein-protein interaction extraction

---

### HGNC (Human Gene Nomenclature Committee)

Official human gene naming authority.

**Initialization:**
```python
from bioservices import HGNC
h = HGNC()
```

**Key Methods:**
- `search(database_or_query=None, query=None, frmt="json")`: search gene symbols/names
- `fetch(database, query, frmt="json")`: retrieve a gene record
- `get_info(frmt="json")`: list the searchable/stored fields

**Use cases:**
- Standardizing human gene names
- Looking up official gene symbols

---

### MyGeneInfo

Gene annotation and query service.

**Initialization:**
```python
from bioservices import MyGeneInfo
m = MyGeneInfo()
```

**Key Methods:**
- `get_genes(ids, ...)` / `get_one_gene(geneid, ...)`: gene annotation by ID
- `get_queries(...)` / `get_one_query(...)`: batch or single free-text query
- `get_metadata()`, `get_taxonomy()`

**Use cases:**
- Batch gene annotation retrieval
- Gene ID conversion

---

## Chemical Compound Resources

### ChEBI (Chemical Entities of Biological Interest)

Dictionary of molecular entities.

**Initialization:**
```python
from bioservices import ChEBI
c = ChEBI()
```

**Key Methods** (REST since bioservices 1.13 — the SOAP interface is gone):
- `getCompleteEntity(chebi_id)`: full entry as a dict-like `ChebiEntity`; accepts
  `"CHEBI:27732"` or `"27732"`
- `getLiteEntity(search, searchCategory="ALL", maximumResults=200, stars="ALL")`: search
- `getCompleteEntityByList(chebi_ids)`: batch retrieval
- `conv(chebi_id, target)`: cross-references for one source, e.g.
  `conv("CHEBI:10102", "KEGG COMPOUND accession")`

`ChebiEntity` exposes `.chebiId`, `.chebiAsciiName`, `.formula`, `.mass`, `.charge`,
`.smiles`, `.inchiKey`, and `.DatabaseLinks` (list of `(accession, source_name)` pairs).
Note `.formula` — the old SOAP attribute was `Formulae`.

**Use cases:**
- Small molecule information
- Chemical structure data
- Compound property lookup

---

### ChEMBL

Bioactive drug-like compound database.

**Initialization:**
```python
from bioservices import ChEMBL
c = ChEMBL()
```

**Key Methods (bioservices >= 1.6.0):**
- `get_molecule(query=None, limit=20, offset=0, filters=None)`: retrieve molecule records
  (by ChEMBL ID, list of IDs, or filters)
- `search_molecule(query)`: Free-text molecule search
- `get_target(query)`: Target information
- `get_similarity(smiles_or_id, similarity)`: Find similar compounds
- `get_substructure(smiles_or_id)`: Substructure search

Note: the pre-1.6.0 helpers `get_compound_by_chemblId()` / `get_molecule_form()`
were removed when the API was simplified in 1.6.0; use `get_molecule()` instead.

**Use cases:**
- Drug discovery data
- Find similar compounds
- Bioactivity information
- Target-compound relationships

---

### UniChem

Chemical identifier mapping service.

**Important (bioservices behavior):** The compound-mapping convenience methods
that previously lived on `bioservices.UniChem` (e.g. `get_compound_id_from_kegg`,
`get_src_compound_ids`, `get_all_compound_ids`) were **dropped from bioservices in
2022** and are no longer available. There is also no KEGG -> ChEMBL route in the
KEGG Web Service. For chemical ID cross-referencing today:

- Use `KEGG.conv("chebi", "compound")` for the supported KEGG -> ChEBI mapping
  (ChEBI IDs are also embedded in KEGG compound entries and can be parsed out).
- For KEGG/ChEBI -> ChEMBL (or other UniChem sources), call the **live UniChem
  REST API directly** (`https://www.ebi.ac.uk/unichem/`) or the **ChEMBL web
  service / `chembl_webresource_client`** — these are separate from bioservices.

**Use cases:**
- Cross-database compound ID mapping (via the routes above, not bioservices.UniChem)
- Linking chemical databases

---

### PubChem

Chemical compound database from NIH.

**Initialization:**
```python
from bioservices import PubChem
p = PubChem()
```

**Key Methods** (PUG REST, refreshed in bioservices 1.14):
- `get_cids_by_name(name)`, `get_cids_by_smiles(smiles)`, `get_cids_by_inchikey(...)`
- `get_compound_by_cid(cid)`, `get_compound_by_name(name)`, `get_compound_by_smiles(...)`
- `get_properties(identifier, namespace="cid", properties=None)`, `get_synonyms(...)`,
  `get_xrefs(...)`, `get_assay(aid)`

**Use cases:**
- Chemical structure retrieval
- Compound property information

---

## Sequence Analysis Tools

### NCBIblast

Sequence similarity searching.

**Initialization:**
```python
from bioservices import NCBIblast
s = NCBIblast(verbose=False)
```

**Key Methods:**
- `run(program, sequence, stype, database, email, **params)`
  - Submit BLAST job
  - `program`: "blastp", "blastn", "blastx", "tblastn", "tblastx"
  - `stype`: "protein" or "dna"
  - `database`: "uniprotkb", "pdb", "refseq_protein", etc.
  - `email`: Required by NCBI
  - Returns: Job ID

- `get_status(jobid)`
  - Check job status
  - Returns: "RUNNING", "FINISHED", "ERROR", "FAILURE", or "NOT_FOUND"

- `wait(jobid)`
  - Block until the job finishes (polls at `checkInterval` seconds)

- `get_result(jobid, result_type)` / `get_result_types(jobid)`
  - Retrieve results; `result_type` is one of the identifiers `get_result_types` returns
    (e.g. "out", "ids", "xml")

> The camelCase `getStatus` / `getResult` / `parametersDetails` names no longer exist;
> only the docstrings still mention them. `NCBIBlastAPI` (bioservices 1.16) offers the same
> run / get_status / get_result flow against NCBI's own BLAST URL API.

**Important:** BLAST jobs are asynchronous. Always check status before retrieving results.

**Use cases:**
- Protein homology searches
- Sequence similarity analysis
- Functional annotation by homology

---

## Pathway & Interaction Resources

### Reactome

Pathway database.

**Initialization:**
```python
from bioservices import Reactome
r = Reactome()
```

**Key Methods:**
- `get_pathways_top(species)`, `get_pathway_containedEvents(identifier)`,
  `get_event_ancestors(identifier)`, `get_complex_subunits(identifier)`
- `search_query(query)`, `search_facet_query(query)`, `get_species_all()`

**Use cases:**
- Human pathway analysis
- Biological process annotation

---

### STRING (replaces PSICQUIC / BioGRID)

`PSICQUIC` and `BioGRID` were **removed in bioservices 1.14**; importing them raises
`ImportError`. STRING covers the same "who interacts with this protein" question.

**Initialization:**
```python
from bioservices import STRING
s = STRING()
```

**Key Methods:**
- `get_interaction_partners(identifiers, species=None, required_score=None, limit=None, network_type="functional")`
  - Partners of the query proteins, including ones outside the input set
- `get_interactions(identifiers, species=...)` — edges *within* the given set
- `get_network(...)`, `get_enrichment(...)`, `get_functional_annotation(...)`,
  `get_ppi_enrichment(...)`, `get_homology(...)`, `get_string_ids(...)`, `get_version()`

**Parameters that matter:** `species` is an NCBI taxid (9606 = human); `required_score`
is 0–1000 (returned `score` values are 0–1); `network_type` is `"functional"` (default)
or `"physical"`.

**Use cases:**
- Protein-protein interaction discovery
- Network analysis and enrichment of an interactor set
- Interactome mapping

---

### IntactComplex

Protein complex database.

**Initialization:**
```python
from bioservices import IntactComplex
i = IntactComplex()
```

**Key Methods:**
- `search(query)`: Search complexes
- `details(complex_ac)`: Complex details

**Use cases:**
- Protein complex composition
- Multi-protein assembly analysis

---

### OmniPath

Integrated signaling pathway database.

**Initialization:**
```python
from bioservices import OmniPath
o = OmniPath()
```

**Key Methods:**
- `get_interactions(query="", frmt="json", fields=[])`
- `get_ptms(query="", ptm_type=None, frmt="json", fields=[])`
- `get_network(frmt="json")`, `get_resources(frmt="json")`

**Use cases:**
- Cell signaling analysis
- Regulatory network mapping

---

## Gene Ontology

### QuickGO

Gene Ontology annotation service.

**Initialization:**
```python
from bioservices import QuickGO
g = QuickGO()
```

**Key Methods** (the QuickGO REST refresh renamed most of these):
- `get_go_terms(query)` / `go_search(query, limit=600, page=1)`
  - Retrieve or search GO term information (parsed JSON)
- `get_go_ancestors(query, relations=...)`, `get_go_children(query)`, `get_go_paths(_from, _to)`
  - Navigate the ontology graph
- `Annotation(geneProductId=None, goId=None, taxonId=None, aspect=None, includeFields=None, limit=100, page=1, ...)`
  - Get GO annotations. `geneProductId` is prefixed (`"UniProtKB:P43403"`), `limit` is
    capped at 100 (higher raises `TypeError`), and the result is a dict with
    `numberOfHits` plus a `results` list of records (`goId`, `goName`, `goAspect`,
    `qualifier`, `evidenceCode`, ...). The old `protein=` / `format=` parameters are gone.
- `Annotation_from_goid(goId, ...)`, `gene_product_search(...)`

**GO categories:**
- Biological Process (BP)
- Molecular Function (MF)
- Cellular Component (CC)

**Use cases:**
- Functional annotation
- Enrichment analysis
- GO term lookup

---

## Genomic Resources

### BioMart

Data mining tool for genomic data.

**Initialization:**
```python
from bioservices import BioMart
b = BioMart()
```

**Key Methods:**
- `registry()`, `datasets(mart)`, `attributes(dataset)`, `filters(dataset)`
- `new_query()` + `add_dataset_to_xml` / `add_attribute_to_xml` / `add_filter_to_xml`
  + `get_xml()`, then `query(xmlq)`

**Use cases:**
- Bulk genomic data retrieval
- Custom genome annotations
- SNP information

---

### ArrayExpress

Gene expression database.

**Initialization:**
```python
from bioservices import ArrayExpress
a = ArrayExpress()
```

**Key Methods:**
- `search(query, page=1, page_size=20, ...)`: search studies (current BioStudies-backed API)
- `get_study(accession)`, `get_files(accession)`, `retrieve_file(accession, filename)`
- The legacy `queryExperiments` / `retrieveExperiment` helpers remain but target the
  retired ArrayExpress endpoints

**Use cases:**
- Gene expression data
- Microarray analysis
- RNA-seq data retrieval

---

### ENA (European Nucleotide Archive)

Nucleotide sequence database.

**Initialization:**
```python
from bioservices import ENA
e = ENA()
```

**Key Methods:**
- `get_data(identifier, frmt=...)`: retrieve records by accession
- `get_taxon(taxon)`, `data_warehouse()`

**Use cases:**
- Nucleotide sequence retrieval
- Genome assembly access

---

## Structural Biology

### PDB (Protein Data Bank)

3D protein structure database.

**Initialization:**
```python
from bioservices import PDB
p = PDB()
```

**Key Methods:**
- `search(query, request_options=None, request_info=None, return_type=None)`: RCSB Search API v2
- `get_current_ids()`, `get_similarity_sequence(seq)`

For downloading coordinate files, fetch from RCSB directly (or use `alterlab-pdb`);
the v2 API wrapper in bioservices is search-oriented.

**Use cases:**
- 3D structure retrieval
- Structure-based analysis
- PyMOL visualization

---

### Pfam

Protein family database.

**Initialization:**
```python
from bioservices import Pfam
p = Pfam()
```

**Key Methods:**
- `show(Id)`, `get_protein(ID, output="json")`

Pfam is now served through InterPro; the bioservices class scrapes those pages rather
than calling a dedicated Pfam REST API, so prefer `alterlab-interpro` for real work.

**Use cases:**
- Protein domain identification
- Family classification
- Functional motif discovery

---

## Specialized Resources

### BioModels

Systems biology model repository.

**Initialization:**
```python
from bioservices import BioModels
b = BioModels()
```

**Key Methods:**
- `get_model(model_id, frmt="json")`, `get_model_files(model_id)`,
  `get_model_download(model_id, filename=...)`, `search(query)`

**Use cases:**
- Systems biology modeling
- SBML model retrieval

---

### COG (Clusters of Orthologous Genes)

Orthologous gene classification.

**Initialization:**
```python
from bioservices import COG
c = COG()
```

**Use cases:**
- Orthology analysis
- Functional classification

---

### BiGG Models

Metabolic network models.

**Initialization:**
```python
from bioservices import BiGG
b = BiGG()
```

**Key Methods:**
- `models` (property): available models
- `get_model(model_id)`, `metabolites(...)`, `reactions(...)`, `genes(model_id)`,
  `search(query, type_)`, `download(model_id, format_="json")`

**Use cases:**
- Metabolic network analysis
- Flux balance analysis

---

## General Patterns

### Error Handling

All services may throw exceptions. Wrap calls in try-except:

```python
try:
    result = service.method(params)
    if result:
        # Process result
        pass
except Exception as e:
    print(f"Error: {e}")
```

### Verbosity Control

Most services support `verbose` parameter:
```python
service = Service(verbose=False)  # Suppress HTTP logs
```

### Rate Limiting

Timeouts live on the `REST` object each service holds (`service.services`), so set them
there — most classes are plain wrappers now rather than `REST` subclasses:
```python
k = KEGG()
k.services.TIMEOUT = 30       # seconds
k.services.settings.TIMEOUT = 30   # equivalent, via the settings object
```

### Output Formats

Common format parameters:
- `frmt`: "xml", "json", "tab", "txt", "fasta"
- `format`: Service-specific variants

### Caching

Caching is opt-in at construction and backed by `requests_cache`:
```python
k = KEGG(cache=True)       # store responses in a local sqlite cache
k.services.clear_cache()   # drop it
```

## Additional Resources

For detailed API documentation:
- Official docs: https://bioservices.readthedocs.io/
- Individual service docs linked from main page
- Source code: https://github.com/cokelaer/bioservices

# BioServices: Common Workflow Patterns

This document describes detailed multi-step workflows for common bioinformatics tasks using BioServices.

## Table of Contents

1. [Complete Protein Analysis Pipeline](#complete-protein-analysis-pipeline)
2. [Pathway Discovery and Network Analysis](#pathway-discovery-and-network-analysis)
3. [Compound Multi-Database Search](#compound-multi-database-search)
4. [Batch Identifier Conversion](#batch-identifier-conversion)
5. [Gene Functional Annotation](#gene-functional-annotation)
6. [Protein Interaction Network Construction](#protein-interaction-network-construction)
7. [Multi-Organism Comparative Analysis](#multi-organism-comparative-analysis)

---

## Complete Protein Analysis Pipeline

**Goal:** Given a protein name, retrieve sequence, find homologs, identify pathways, and discover interactions.

**Example:** Analyzing human ZAP70 protein

> Two API facts these workflows depend on (bioservices 1.16):
> `UniProt.search(frmt=...)` accepts only `tsv`/`xlsx`/`fasta`/`json`/`gff` with
> UniProt return-field names in `columns`, and `UniProt.mapping()` returns
> `{"results": [{"from", "to"}...], "failedIds": [...]}`. The snippets below use this
> helper to collapse a mapping job into `{source_id: [target_ids]}`:
>
> ```python
> from collections import defaultdict
>
> def as_dict(job):
>     out = defaultdict(list)
>     for row in (job or {}).get("results", []):
>         out[row["from"]].append(row["to"])
>     return dict(out)
> ```

### Step 1: UniProt Search and Identifier Retrieval

```python
from bioservices import UniProt

u = UniProt(verbose=False)

# Search for protein by name
query = "ZAP70_HUMAN"
results = u.search(query, frmt="tsv", columns="accession,gene_names,organism_name,length")

# Parse results
lines = results.strip().split("\n")
if len(lines) > 1:
    header = lines[0]
    data = lines[1].split("\t")
    uniprot_id = data[0]  # e.g., P43403
    gene_names = data[1]   # e.g., ZAP70

print(f"UniProt ID: {uniprot_id}")
print(f"Gene names: {gene_names}")
```

**Output:**
- UniProt accession: P43403
- Gene name: ZAP70

### Step 2: Sequence Retrieval

```python
# Retrieve FASTA sequence
sequence = u.retrieve(uniprot_id, frmt="fasta")
print(sequence)

# Extract just the sequence string (remove header)
seq_lines = sequence.split("\n")
sequence_only = "".join(seq_lines[1:])  # Skip FASTA header
```

**Output:** Complete protein sequence in FASTA format

### Step 3: BLAST Similarity Search

```python
from bioservices import NCBIblast
import time

s = NCBIblast(verbose=False)

# Submit BLAST job
jobid = s.run(
    program="blastp",
    sequence=sequence_only,
    stype="protein",
    database="uniprotkb",
    email="your.email@example.com"
)

print(f"BLAST Job ID: {jobid}")

# Wait for completion (method names are snake_case in current bioservices)
while True:
    status = s.get_status(jobid)
    print(f"Status: {status}")
    if status == "FINISHED":
        break
    elif status in ("ERROR", "FAILURE", "NOT_FOUND"):
        print(f"BLAST job did not complete: {status}")
        break
    time.sleep(5)

# Retrieve results
if status == "FINISHED":
    blast_results = s.get_result(jobid, "out")
    print(blast_results[:500])  # Print first 500 characters
```

**Output:** BLAST alignment results showing similar proteins

### Step 4: KEGG Pathway Discovery

```python
from bioservices import KEGG

k = KEGG()

# Get KEGG gene ID from UniProt mapping
kegg_mapping = as_dict(u.mapping(fr="UniProtKB_AC-ID", to="KEGG", query=uniprot_id))
print(f"KEGG mapping: {kegg_mapping}")

# Extract KEGG gene ID (e.g., hsa:7535)
if kegg_mapping:
    kegg_gene_id = kegg_mapping.get(uniprot_id, [None])[0]

    if kegg_gene_id:
        # Find pathways containing this gene
        organism = kegg_gene_id.split(":")[0]  # e.g., "hsa"
        gene_id = kegg_gene_id.split(":")[1]   # e.g., "7535"

        # get_pathway_by_gene returns {pathway_id: name} — the names come free.
        pathways = k.get_pathway_by_gene(gene_id, organism)
        print(f"Found {len(pathways)} pathways:")

        for pathway_id, pathway_name in pathways.items():
            print(f"  {pathway_id}: {pathway_name}")
```

**Output:**
- hsa04064 - NF-kappa B signaling pathway
- hsa04650 - Natural killer cell mediated cytotoxicity
- hsa04660 - T cell receptor signaling pathway
- hsa04662 - B cell receptor signaling pathway

### Step 5: Protein-Protein Interactions

```python
from bioservices import STRING   # PSICQUIC/BioGRID were removed in bioservices 1.14

s = STRING()

# Human (taxid 9606) partners of ZAP70, medium confidence and above
interactions = s.get_interaction_partners("ZAP70", species=9606, required_score=400, limit=25)

print(f"Found {len(interactions)} interactions")
for row in interactions[:5]:
    print(f"  {row['preferredName_A']} - {row['preferredName_B']}: score={row['score']}")
```

**Output:** List of proteins that interact with ZAP70

### Step 6: Gene Ontology Annotation

```python
from bioservices import QuickGO

g = QuickGO()

# Get GO annotations for protein. The QuickGO REST parameters are
# geneProductId / goId / taxonId / aspect; limit is capped at 100 per page.
annotations = g.Annotation(
    geneProductId=f"UniProtKB:{uniprot_id}",
    includeFields="goName",
    limit=100,
)

if annotations:
    print(f"Found {annotations['numberOfHits']} GO annotations")
    for row in annotations["results"][:5]:
        print(f"  {row['goId']}: {row['goName']} [{row['goAspect']}]")
```

**Output:** GO terms annotating ZAP70 function, process, and location

### Complete Pipeline Summary

**Inputs:** Protein name (e.g., "ZAP70_HUMAN")

**Outputs:**
1. UniProt accession and gene name
2. Protein sequence (FASTA)
3. Similar proteins (BLAST results)
4. Biological pathways (KEGG)
5. Interaction partners (STRING)
6. Functional annotations (GO terms)

**Script:** `scripts/protein_analysis_workflow.py` automates this entire pipeline.

---

## Pathway Discovery and Network Analysis

**Goal:** Analyze all pathways for an organism and extract protein interaction networks.

**Example:** Human (hsa) pathway analysis

### Step 1: Get All Pathways for Organism

```python
from bioservices import KEGG

k = KEGG()
k.organism = "hsa"

# Get all pathway IDs
pathway_ids = k.pathwayIds
print(f"Found {len(pathway_ids)} pathways for {k.organism}")

# Display first few
for pid in pathway_ids[:10]:
    print(f"  {pid}")
```

**Output:** List of ~300 human pathways

### Step 2: Parse Pathway for Interactions

```python
# Analyze specific pathway
pathway_id = "hsa04660"  # T cell receptor signaling

# Get KGML data
kgml_data = k.parse_kgml_pathway(pathway_id)

# Extract entries (genes/proteins)
entries = kgml_data['entries']
print(f"Pathway contains {len(entries)} entries")

# Extract relations (interactions)
relations = kgml_data['relations']
print(f"Found {len(relations)} relations")

# Analyze relation types
relation_types = {}
for rel in relations:
    rel_type = rel.get('name', 'unknown')
    relation_types[rel_type] = relation_types.get(rel_type, 0) + 1

print("\nRelation type distribution:")
for rel_type, count in sorted(relation_types.items()):
    print(f"  {rel_type}: {count}")
```

**Output:**
- Entry count (genes/proteins in pathway)
- Relation count (interactions)
- Distribution of interaction types (activation, inhibition, binding, etc.)

### Step 3: Extract Protein-Protein Interactions

```python
# Filter for specific interaction types
pprel_interactions = [
    rel for rel in relations
    if rel.get('link') == 'PPrel'  # Protein-protein relation
]

print(f"Found {len(pprel_interactions)} protein-protein interactions")

# Extract interaction details
for rel in pprel_interactions[:10]:
    entry1 = rel['entry1']
    entry2 = rel['entry2']
    interaction_type = rel.get('name', 'unknown')

    print(f"  {entry1} -> {entry2}: {interaction_type}")
```

**Output:** Directed protein-protein interactions with types

### Step 4: Convert to Network Format (SIF)

```python
# Get Simple Interaction Format (filters for key interactions)
sif_data = k.pathway2sif(pathway_id)

# SIF format: source, interaction_type, target
print("\nSimple Interaction Format:")
for interaction in sif_data[:10]:
    print(f"  {interaction}")
```

**Output:** Network edges suitable for Cytoscape or NetworkX

### Step 5: Batch Analysis of All Pathways

```python
import pandas as pd

# Analyze all pathways (this takes time!)
all_results = []

for pathway_id in pathway_ids[:50]:  # Limit for example
    try:
        kgml = k.parse_kgml_pathway(pathway_id)

        result = {
            'pathway_id': pathway_id,
            'num_entries': len(kgml.get('entries', [])),
            'num_relations': len(kgml.get('relations', []))
        }

        all_results.append(result)

    except Exception as e:
        print(f"Error parsing {pathway_id}: {e}")

# Create DataFrame
df = pd.DataFrame(all_results)
print(df.describe())

# Find largest pathways
print("\nLargest pathways:")
print(df.nlargest(10, 'num_entries')[['pathway_id', 'num_entries', 'num_relations']])
```

**Output:** Statistical summary of pathway sizes and interaction densities

**Script:** `scripts/pathway_analysis.py` implements this workflow with export options.

---

## Compound Multi-Database Search

**Goal:** Search for compound by name and retrieve identifiers across KEGG, ChEBI, and ChEMBL.

**Example:** Geldanamycin (antibiotic)

### Step 1: Search KEGG Compound Database

```python
from bioservices import KEGG

k = KEGG()

# Search by compound name
compound_name = "Geldanamycin"
results = k.find("compound", compound_name)

print(f"KEGG search results for '{compound_name}':")
print(results)

# Extract compound ID
if results:
    lines = results.strip().split("\n")
    if lines:
        kegg_id = lines[0].split("\t")[0]  # e.g., cpd:C11222
        kegg_id_clean = kegg_id.replace("cpd:", "")  # C11222
        print(f"\nKEGG Compound ID: {kegg_id_clean}")
```

**Output:** KEGG ID (e.g., C11222)

### Step 2: Get KEGG Entry with Database Links

```python
# Retrieve compound entry
compound_entry = k.get(kegg_id)

# Parse entry for database links
chebi_id = None
for line in compound_entry.split("\n"):
    if "ChEBI:" in line:
        # Extract ChEBI ID
        parts = line.split("ChEBI:")
        if len(parts) > 1:
            chebi_id = parts[1].strip().split()[0]
            print(f"ChEBI ID: {chebi_id}")
            break

# Display entry snippet
print("\nKEGG Entry (first 500 chars):")
print(compound_entry[:500])
```

**Output:** ChEBI ID (e.g., 5292) and compound information

### Step 3: Confirm KEGG -> ChEBI via KEGG.conv

```python
# KEGG.conv is the supported bioservices route for KEGG -> ChEBI.
# It returns a dict keyed by source IDs, e.g. {'cpd:C11222': 'chebi:5292', ...}.
conv = k.conv("chebi", "compound")
chebi_from_conv = conv.get(kegg_id)  # kegg_id is the 'cpd:Cxxxxx' form
print(f"ChEBI (via conv): {chebi_from_conv}")
```

> **No bioservices KEGG -> ChEMBL route.** The old `UniChem.get_compound_id_from_kegg`
> helper was dropped from bioservices in 2022, and the KEGG Web Service has no
> KEGG -> ChEMBL mapping. If you need a ChEMBL ID, fetch it separately via the
> ChEMBL web service / `chembl_webresource_client` or the live UniChem REST API
> (`https://www.ebi.ac.uk/unichem/`) — both are outside bioservices.

### Step 4: Retrieve Detailed ChEBI Information

```python
# Get ChEBI information for the cross-referenced compound
if chebi_id:
    from bioservices import ChEBI
    c = ChEBI()

    try:
        # REST-backed since bioservices 1.13: the entity is dict-like, and the
        # formula attribute is `.formula` (SOAP-era code used `.Formulae`).
        chebi_entity = c.getCompleteEntity(f"CHEBI:{chebi_id}")
        print(f"\nChEBI Formula: {chebi_entity.formula}")
        print(f"ChEBI Name: {chebi_entity.chebiAsciiName}")
        print(f"ChEBI SMILES: {chebi_entity.smiles}")
    except Exception as e:
        print(f"ChEBI lookup failed: {e}")
```

**Output:** Chemical properties from KEGG and ChEBI

### Complete Compound Workflow Summary

**Input:** Compound name (e.g., "Geldanamycin")

**Output:**
- KEGG ID: C11222
- ChEBI ID: 5292 (via KEGG entry / `KEGG.conv`)
- Chemical formula
- Molecular weight
- ChEMBL ID / SMILES: not via bioservices — fetch from ChEMBL separately if needed

**Script:** `scripts/compound_cross_reference.py` automates the KEGG -> ChEBI portion.

---

## Batch Identifier Conversion

**Goal:** Convert multiple identifiers between databases efficiently.

### Batch UniProt → KEGG Mapping

```python
from bioservices import UniProt

u = UniProt()

# List of UniProt IDs
uniprot_ids = ["P43403", "P04637", "P53779", "Q9Y6K9"]

# Batch mapping (comma-separated)
query_string = ",".join(uniprot_ids)
results = as_dict(u.mapping(fr="UniProtKB_AC-ID", to="KEGG", query=query_string))

print("UniProt → KEGG mapping:")
for uniprot_id, kegg_ids in results.items():
    print(f"  {uniprot_id} → {kegg_ids}")
```

**Output:** Dictionary mapping each UniProt ID to KEGG gene IDs

### Batch File Processing

```python
import csv

# Read identifiers from file
def read_ids_from_file(filename):
    with open(filename, 'r') as f:
        ids = [line.strip() for line in f if line.strip()]
    return ids

# Process in chunks (API limits)
def batch_convert(ids, from_db, to_db, chunk_size=100):
    u = UniProt()
    all_results = {}

    for i in range(0, len(ids), chunk_size):
        chunk = ids[i:i+chunk_size]
        query = ",".join(chunk)

        try:
            all_results.update(as_dict(u.mapping(fr=from_db, to=to_db, query=query)))
            print(f"Processed {min(i+chunk_size, len(ids))}/{len(ids)}")
        except Exception as e:
            print(f"Error processing chunk {i}: {e}")

    return all_results

# Write results to CSV
def write_mapping_to_csv(mapping, output_file):
    with open(output_file, 'w', newline='') as f:
        writer = csv.writer(f)
        writer.writerow(['Source_ID', 'Target_IDs'])

        for source_id, target_ids in mapping.items():
            target_str = ";".join(target_ids) if target_ids else "No mapping"
            writer.writerow([source_id, target_str])

# Example usage
input_ids = read_ids_from_file("uniprot_ids.txt")
mapping = batch_convert(input_ids, "UniProtKB_AC-ID", "KEGG", chunk_size=50)
write_mapping_to_csv(mapping, "uniprot_to_kegg_mapping.csv")
```

**Script:** `scripts/batch_id_converter.py` provides command-line batch conversion.

---

## Gene Functional Annotation

**Goal:** Retrieve comprehensive functional information for a gene.

### Workflow

```python
from bioservices import UniProt, KEGG, QuickGO

# Gene of interest
gene_symbol = "TP53"

# 1. Find UniProt entry
u = UniProt()
search_results = u.search(f"gene:{gene_symbol} AND organism_id:9606",
                          frmt="tsv",
                          columns="accession,gene_names,protein_name")

# Extract UniProt ID
lines = search_results.strip().split("\n")
if len(lines) > 1:
    uniprot_id = lines[1].split("\t")[0]
    protein_name = lines[1].split("\t")[2]
    print(f"Protein: {protein_name}")
    print(f"UniProt ID: {uniprot_id}")

# 2. Get KEGG pathways
kegg_mapping = as_dict(u.mapping(fr="UniProtKB_AC-ID", to="KEGG", query=uniprot_id))
if uniprot_id in kegg_mapping:
    kegg_id = kegg_mapping[uniprot_id][0]

    k = KEGG()
    organism, gene_id = kegg_id.split(":")
    pathways = k.get_pathway_by_gene(gene_id, organism)

    print(f"\nPathways ({len(pathways)}):")
    for pathway_id in pathways[:5]:
        print(f"  {pathway_id}")

# 3. Get GO annotations
g = QuickGO()
go_annotations = g.Annotation(geneProductId=f"UniProtKB:{uniprot_id}", includeFields="goName", limit=100)

if go_annotations:
    print(f"\nGO Annotations ({go_annotations['numberOfHits']} total):")

    # Group by aspect; goAspect is spelled out (biological_process, ...)
    aspects = {}
    for row in go_annotations["results"]:
        aspects.setdefault(row["goAspect"], []).append(row.get("goName", row["goId"]))

    for aspect, terms in aspects.items():
        print(f"  {aspect}: {len(terms)} terms")

# 4. Get protein sequence features
full_entry = u.retrieve(uniprot_id, frmt="txt")
print("\nProtein Features:")
for line in full_entry.split("\n"):
    if line.startswith("FT   DOMAIN"):
        print(f"  {line}")
```

**Output:** Comprehensive annotation including name, pathways, GO terms, and features.

---

## Protein Interaction Network Construction

**Goal:** Build a protein-protein interaction network for a set of proteins.

### Workflow

```python
from bioservices import STRING
import networkx as nx

# Proteins of interest
proteins = ["ZAP70", "LCK", "LAT", "LCP2", "PLCG1"]

s = STRING()
G = nx.Graph()

# One call covers the whole set; STRING returns the edges among them.
try:
    for row in s.get_interactions(proteins, species=9606, required_score=400):
        G.add_edge(row["preferredName_A"], row["preferredName_B"], weight=row["score"])
except Exception as e:
    print(f"STRING query failed: {e}")

# To pull in partners outside the input set, add:
#   s.get_interaction_partners(proteins, species=9606, limit=10)

print(f"Network: {G.number_of_nodes()} nodes, {G.number_of_edges()} edges")

# Analyze network
print("\nNode degrees:")
for node in proteins:
    if node in G:
        print(f"  {node}: {G.degree(node)} interactions")

# Export for visualization
nx.write_gml(G, "protein_network.gml")
print("\nNetwork exported to protein_network.gml")
```

**Output:** NetworkX graph exported in GML format for Cytoscape visualization.

---

## Multi-Organism Comparative Analysis

**Goal:** Compare pathway or gene presence across multiple organisms.

### Workflow

```python
from bioservices import KEGG

k = KEGG()

# Organisms to compare
organisms = ["hsa", "mmu", "dme", "sce"]  # Human, mouse, fly, yeast
organism_names = {
    "hsa": "Human",
    "mmu": "Mouse",
    "dme": "Fly",
    "sce": "Yeast"
}

# Pathway of interest
pathway_name = "cell cycle"

print(f"Searching for '{pathway_name}' pathway across organisms:\n")

for org in organisms:
    k.organism = org

    # Search pathways
    results = k.lookfor_pathway(pathway_name)

    print(f"{organism_names[org]} ({org}):")
    if results:
        for pathway in results[:3]:  # Show first 3
            print(f"  {pathway}")
    else:
        print("  No matches found")
    print()
```

**Output:** Pathway presence/absence across organisms.

---

## Best Practices for Workflows

### 1. Error Handling

Always wrap service calls:
```python
try:
    result = service.method(params)
    if result:
        # Process
        pass
except Exception as e:
    print(f"Error: {e}")
```

### 2. Rate Limiting

Add delays for batch processing:
```python
import time

for item in items:
    result = service.query(item)
    time.sleep(0.5)  # 500ms delay
```

### 3. Result Validation

Check for empty or unexpected results:
```python
if result and len(result) > 0:
    # Process
    pass
else:
    print("No results returned")
```

### 4. Progress Reporting

For long workflows:
```python
total = len(items)
for i, item in enumerate(items):
    # Process item
    if (i + 1) % 10 == 0:
        print(f"Processed {i+1}/{total}")
```

### 5. Data Export

Save intermediate results:
```python
import json

with open("results.json", "w") as f:
    json.dump(results, f, indent=2)
```

---

## Integration with Other Tools

### BioPython Integration

```python
from bioservices import UniProt
from Bio import SeqIO
from io import StringIO

u = UniProt()
fasta_data = u.retrieve("P43403", "fasta")

# Parse with BioPython
fasta_io = StringIO(fasta_data)
record = SeqIO.read(fasta_io, "fasta")

print(f"Sequence length: {len(record.seq)}")
print(f"Description: {record.description}")
```

### Pandas Integration

```python
from bioservices import UniProt
import pandas as pd
from io import StringIO

u = UniProt()
results = u.search("zap70", frmt="tsv", columns="accession,gene_names,length,organism_name")

# Load into DataFrame
df = pd.read_csv(StringIO(results), sep="\t")
print(df.head())
print(df.describe())
```

### NetworkX Integration

See Protein Interaction Network Construction above.

---

For complete working examples, see the scripts in `scripts/` directory.

#!/usr/bin/env python3
"""
Batch Identifier Converter

This script converts multiple identifiers between biological databases
using UniProt's mapping service. Supports batch processing with
automatic chunking and error handling.

Usage:
    python batch_id_converter.py INPUT_FILE --from DB1 --to DB2 [options]

Examples:
    python batch_id_converter.py uniprot_ids.txt --from UniProtKB_AC-ID --to KEGG
    python batch_id_converter.py gene_ids.txt --from GeneID --to UniProtKB --output mapping.csv
    python batch_id_converter.py ids.txt --from UniProtKB_AC-ID --to Ensembl --chunk-size 50

Input file format:
    One identifier per line (plain text)

Common database codes:
    UniProtKB_AC-ID  - UniProt accession/ID
    KEGG             - KEGG gene IDs
    GeneID           - NCBI Gene (Entrez) IDs
    Ensembl          - Ensembl gene IDs
    Ensembl_Protein  - Ensembl protein IDs
    RefSeq_Protein   - RefSeq protein IDs
    PDB              - Protein Data Bank IDs
    HGNC             - Human gene symbols
    GO               - Gene Ontology IDs
"""

import sys
import argparse
import csv
import time
from bioservices import UniProt


# Common database code mappings
DATABASE_CODES = {
    'uniprot': 'UniProtKB_AC-ID',
    'uniprotkb': 'UniProtKB_AC-ID',
    'kegg': 'KEGG',
    'geneid': 'GeneID',
    'entrez': 'GeneID',
    'ensembl': 'Ensembl',
    'ensembl_protein': 'Ensembl_Protein',
    'ensembl_transcript': 'Ensembl_Transcript',
    'refseq': 'RefSeq_Protein',
    'refseq_protein': 'RefSeq_Protein',
    'pdb': 'PDB',
    'hgnc': 'HGNC',
    'mgi': 'MGI',
    'go': 'GO',
    'pfam': 'Pfam',
    'interpro': 'InterPro',
    'reactome': 'Reactome',
    'string': 'STRING',
    'biogrid': 'BioGRID'
}

# Note: the legal (from, to) pairs are published by UniProt itself — inspect
# UniProt().valid_mapping if a pair is rejected. Non-UniProt sources such as
# KEGG or PDB map only *to* UniProtKB / UniProtKB-Swiss-Prot / UniParc.


def _mapping_as_dict(job):
    """Collapse a UniProt ID-mapping job into {source_id: [target_ids]}.

    Since the June-2022 UniProt API (bioservices >= 1.10), ``UniProt.mapping()``
    returns the raw job payload — ``{"results": [{"from": ..., "to": ...}, ...],
    "failedIds": [...]}`` — rather than a dict keyed by the source identifier.
    Some targets (e.g. ``to="UniProtKB"``) return an entry object instead of a
    plain accession, so those are reduced to their ``primaryAccession``.
    """
    out = {}
    for row in (job or {}).get("results", []):
        target = row.get("to")
        if isinstance(target, dict):
            target = target.get("primaryAccession") or target.get("id") or str(target)
        out.setdefault(row.get("from"), []).append(target)
    return out


def normalize_database_code(code):
    """Normalize database code to official format."""
    # Try exact match first
    if code in DATABASE_CODES.values():
        return code

    # Try lowercase lookup
    lowercase = code.lower()
    if lowercase in DATABASE_CODES:
        return DATABASE_CODES[lowercase]

    # Return as-is if not found (may still be valid)
    return code


def read_ids_from_file(filename):
    """Read identifiers from file (one per line)."""
    print(f"Reading identifiers from {filename}...")

    ids = []
    with open(filename, 'r') as f:
        for line in f:
            line = line.strip()
            if line and not line.startswith('#'):
                ids.append(line)

    print(f"✓ Read {len(ids)} identifier(s)")

    return ids


def batch_convert(ids, from_db, to_db, chunk_size=100, delay=0.5):
    """Convert IDs with automatic chunking and error handling."""
    print(f"\nConverting {len(ids)} IDs:")
    print(f"  From: {from_db}")
    print(f"  To: {to_db}")
    print(f"  Chunk size: {chunk_size}")
    print()

    u = UniProt(verbose=False)
    all_results = {}
    failed_ids = []

    total_chunks = (len(ids) + chunk_size - 1) // chunk_size

    for i in range(0, len(ids), chunk_size):
        chunk = ids[i:i+chunk_size]
        chunk_num = (i // chunk_size) + 1

        query = ",".join(chunk)

        try:
            print(f"  [{chunk_num}/{total_chunks}] Processing {len(chunk)} IDs...", end=" ")

            results = _mapping_as_dict(u.mapping(fr=from_db, to=to_db, query=query))

            if results:
                all_results.update(results)
                mapped_count = len([v for v in results.values() if v])
                print(f"✓ Mapped: {mapped_count}/{len(chunk)}")
            else:
                print(f"✗ No mappings returned")
                failed_ids.extend(chunk)

            # Rate limiting
            if delay > 0 and i + chunk_size < len(ids):
                time.sleep(delay)

        except Exception as e:
            print(f"✗ Error: {e}")

            # Try individual IDs in failed chunk
            print(f"    Retrying individual IDs...")
            for single_id in chunk:
                try:
                    result = _mapping_as_dict(u.mapping(fr=from_db, to=to_db, query=single_id))
                    if result:
                        all_results.update(result)
                        print(f"      ✓ {single_id}")
                    else:
                        failed_ids.append(single_id)
                        print(f"      ✗ {single_id} - no mapping")
                except Exception as e2:
                    failed_ids.append(single_id)
                    print(f"      ✗ {single_id} - {e2}")

                time.sleep(0.2)

    # Add missing IDs to results (mark as failed)
    for id_ in ids:
        if id_ not in all_results:
            all_results[id_] = None

    print(f"\n✓ Conversion complete:")
    print(f"  Total: {len(ids)}")
    print(f"  Mapped: {len([v for v in all_results.values() if v])}")
    print(f"  Failed: {len(failed_ids)}")

    return all_results, failed_ids


def save_mapping_csv(mapping, output_file, from_db, to_db):
    """Save mapping results to CSV."""
    print(f"\nSaving results to {output_file}...")

    with open(output_file, 'w', newline='') as f:
        writer = csv.writer(f)

        # Header
        writer.writerow(['Source_ID', 'Source_DB', 'Target_IDs', 'Target_DB', 'Mapping_Status'])

        # Data
        for source_id, target_ids in sorted(mapping.items()):
            if target_ids:
                target_str = ";".join(target_ids)
                status = "Success"
            else:
                target_str = ""
                status = "Failed"

            writer.writerow([source_id, from_db, target_str, to_db, status])

    print(f"✓ Results saved")


def save_failed_ids(failed_ids, output_file):
    """Save failed IDs to file."""
    if not failed_ids:
        return

    print(f"\nSaving failed IDs to {output_file}...")

    with open(output_file, 'w') as f:
        for id_ in failed_ids:
            f.write(f"{id_}\n")

    print(f"✓ Saved {len(failed_ids)} failed ID(s)")


def print_mapping_summary(mapping, from_db, to_db):
    """Print summary of mapping results."""
    print(f"\n{'='*70}")
    print("MAPPING SUMMARY")
    print(f"{'='*70}")

    total = len(mapping)
    mapped = len([v for v in mapping.values() if v])
    failed = total - mapped

    print(f"\nSource database: {from_db}")
    print(f"Target database: {to_db}")
    print(f"\nTotal identifiers: {total}")
    print(f"Successfully mapped: {mapped} ({mapped/total*100:.1f}%)")
    print(f"Failed to map: {failed} ({failed/total*100:.1f}%)")

    # Show some examples
    if mapped > 0:
        print(f"\nExample mappings (first 5):")
        count = 0
        for source_id, target_ids in mapping.items():
            if target_ids:
                target_str = ", ".join(target_ids[:3])
                if len(target_ids) > 3:
                    target_str += f" ... +{len(target_ids)-3} more"
                print(f"  {source_id} → {target_str}")
                count += 1
                if count >= 5:
                    break

    # Show multiple mapping statistics
    multiple_mappings = [v for v in mapping.values() if v and len(v) > 1]
    if multiple_mappings:
        print(f"\nMultiple target mappings: {len(multiple_mappings)} ID(s)")
        print(f"  (These source IDs map to multiple target IDs)")

    print(f"{'='*70}")


def list_common_databases():
    """Print list of common database codes."""
    print("\nCommon Database Codes:")
    print("-" * 70)
    print(f"{'Alias':<20} {'Official Code':<30}")
    print("-" * 70)

    for alias, code in sorted(DATABASE_CODES.items()):
        if alias != code.lower():
            print(f"{alias:<20} {code:<30}")

    print("-" * 70)
    print("\nNote: Many other database codes are supported.")
    print("See UniProt documentation for complete list.")


def main():
    """Main conversion workflow."""
    parser = argparse.ArgumentParser(
        description="Batch convert biological identifiers between databases",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog="""
Examples:
  python batch_id_converter.py uniprot_ids.txt --from UniProtKB_AC-ID --to KEGG
  python batch_id_converter.py ids.txt --from GeneID --to UniProtKB -o mapping.csv
  python batch_id_converter.py ids.txt --from uniprot --to ensembl --chunk-size 50

Common database codes:
  UniProtKB_AC-ID, KEGG, GeneID, Ensembl, Ensembl_Protein,
  RefSeq_Protein, PDB, HGNC, GO, Pfam, InterPro, Reactome

Use --list-databases to see all supported aliases.
        """
    )
    parser.add_argument("input_file", help="Input file with IDs (one per line)")
    parser.add_argument("--from", dest="from_db", required=True,
                       help="Source database code")
    parser.add_argument("--to", dest="to_db", required=True,
                       help="Target database code")
    parser.add_argument("-o", "--output", default=None,
                       help="Output CSV file (default: mapping_results.csv)")
    parser.add_argument("--chunk-size", type=int, default=100,
                       help="Number of IDs per batch (default: 100)")
    parser.add_argument("--delay", type=float, default=0.5,
                       help="Delay between batches in seconds (default: 0.5)")
    parser.add_argument("--save-failed", action="store_true",
                       help="Save failed IDs to separate file")
    parser.add_argument("--list-databases", action="store_true",
                       help="List common database codes and exit")

    args = parser.parse_args()

    # List databases and exit
    if args.list_databases:
        list_common_databases()
        sys.exit(0)

    print("=" * 70)
    print("BIOSERVICES: Batch Identifier Converter")
    print("=" * 70)

    # Normalize database codes
    from_db = normalize_database_code(args.from_db)
    to_db = normalize_database_code(args.to_db)

    if from_db != args.from_db:
        print(f"\nNote: Normalized '{args.from_db}' → '{from_db}'")
    if to_db != args.to_db:
        print(f"Note: Normalized '{args.to_db}' → '{to_db}'")

    # Read input IDs
    try:
        ids = read_ids_from_file(args.input_file)
    except Exception as e:
        print(f"\n✗ Error reading input file: {e}")
        sys.exit(1)

    if not ids:
        print("\n✗ No IDs found in input file")
        sys.exit(1)

    # Perform conversion
    mapping, failed_ids = batch_convert(
        ids,
        from_db,
        to_db,
        chunk_size=args.chunk_size,
        delay=args.delay
    )

    # Print summary
    print_mapping_summary(mapping, from_db, to_db)

    # Save results
    output_file = args.output or "mapping_results.csv"
    save_mapping_csv(mapping, output_file, from_db, to_db)

    # Save failed IDs if requested
    if args.save_failed and failed_ids:
        failed_file = output_file.replace(".csv", "_failed.txt")
        save_failed_ids(failed_ids, failed_file)

    print(f"\n✓ Done!")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Compound Cross-Database Search

This script searches for a compound by name and retrieves identifiers
from multiple databases:
- KEGG Compound
- ChEBI (cross-referenced from the KEGG entry)
- Basic compound properties

Note: bioservices dropped its UniChem KEGG->ChEMBL mapping helper in 2022,
and there is no KEGG Web Service KEGG->ChEMBL route either. If a ChEMBL ID
is required, obtain it separately (e.g. the ChEMBL REST API /
chembl_webresource_client, or the live UniChem REST API directly). This
script therefore stops at the reliable KEGG -> ChEBI cross-reference.

Usage:
    python compound_cross_reference.py COMPOUND_NAME [--output FILE]

Examples:
    python compound_cross_reference.py Geldanamycin
    python compound_cross_reference.py "Adenosine triphosphate"
    python compound_cross_reference.py Aspirin --output aspirin_info.txt
"""

import sys
import argparse
from bioservices import KEGG, ChEBI


def search_kegg_compound(compound_name):
    """Search KEGG for compound by name."""
    print(f"\n{'='*70}")
    print("STEP 1: KEGG Compound Search")
    print(f"{'='*70}")

    k = KEGG()

    print(f"Searching KEGG for: {compound_name}")

    try:
        results = k.find("compound", compound_name)

        if not results or not results.strip():
            print(f"✗ No results found in KEGG")
            return k, None

        # Parse results
        lines = results.strip().split("\n")
        print(f"✓ Found {len(lines)} result(s):\n")

        for i, line in enumerate(lines[:5], 1):
            parts = line.split("\t")
            kegg_id = parts[0]
            description = parts[1] if len(parts) > 1 else "No description"
            print(f"  {i}. {kegg_id}: {description}")

        # Use first result
        first_result = lines[0].split("\t")
        kegg_id = first_result[0].replace("cpd:", "")

        print(f"\nUsing: {kegg_id}")

        return k, kegg_id

    except Exception as e:
        print(f"✗ Error: {e}")
        return k, None


def get_kegg_info(kegg, kegg_id):
    """Retrieve detailed KEGG compound information."""
    print(f"\n{'='*70}")
    print("STEP 2: KEGG Compound Details")
    print(f"{'='*70}")

    try:
        print(f"Retrieving KEGG entry for {kegg_id}...")

        entry = kegg.get(f"cpd:{kegg_id}")

        if not entry:
            print("✗ Failed to retrieve entry")
            return None

        # Parse entry
        compound_info = {
            'kegg_id': kegg_id,
            'name': None,
            'formula': None,
            'exact_mass': None,
            'mol_weight': None,
            'chebi_id': None,
            'pathways': []
        }

        current_section = None

        for line in entry.split("\n"):
            if line.startswith("NAME"):
                compound_info['name'] = line.replace("NAME", "").strip().rstrip(";")

            elif line.startswith("FORMULA"):
                compound_info['formula'] = line.replace("FORMULA", "").strip()

            elif line.startswith("EXACT_MASS"):
                compound_info['exact_mass'] = line.replace("EXACT_MASS", "").strip()

            elif line.startswith("MOL_WEIGHT"):
                compound_info['mol_weight'] = line.replace("MOL_WEIGHT", "").strip()

            elif "ChEBI:" in line:
                parts = line.split("ChEBI:")
                if len(parts) > 1:
                    compound_info['chebi_id'] = parts[1].strip().split()[0]

            elif line.startswith("PATHWAY"):
                current_section = "pathway"
                pathway = line.replace("PATHWAY", "").strip()
                if pathway:
                    compound_info['pathways'].append(pathway)

            elif current_section == "pathway" and line.startswith("            "):
                pathway = line.strip()
                if pathway:
                    compound_info['pathways'].append(pathway)

            elif line.startswith(" ") and not line.startswith("            "):
                current_section = None

        # Display information
        print(f"\n✓ KEGG Compound Information:")
        print(f"  ID: {compound_info['kegg_id']}")
        print(f"  Name: {compound_info['name']}")
        print(f"  Formula: {compound_info['formula']}")
        print(f"  Exact Mass: {compound_info['exact_mass']}")
        print(f"  Molecular Weight: {compound_info['mol_weight']}")

        if compound_info['chebi_id']:
            print(f"  ChEBI ID: {compound_info['chebi_id']}")

        if compound_info['pathways']:
            print(f"  Pathways: {len(compound_info['pathways'])} found")

        return compound_info

    except Exception as e:
        print(f"✗ Error: {e}")
        return None


def confirm_chebi_via_conv(kegg, kegg_id, chebi_from_entry):
    """Confirm/obtain the ChEBI ID via KEGG.conv as a cross-check.

    KEGG.conv('chebi', 'compound') returns a dict keyed by 'cpd:Cxxxxx'.
    This is the supported bioservices route for KEGG -> ChEBI.
    """
    print(f"\n{'='*70}")
    print("STEP 3: KEGG -> ChEBI (via KEGG.conv)")
    print(f"{'='*70}")

    try:
        print(f"Mapping cpd:{kegg_id} to ChEBI via KEGG.conv...")
        conv = kegg.conv("chebi", "compound")
        chebi = conv.get(f"cpd:{kegg_id}")

        if chebi:
            print(f"✓ ChEBI (conv): {chebi}")
            return chebi.replace("chebi:", "")
        else:
            print("⊘ No conv mapping; falling back to ChEBI parsed from entry")
            return chebi_from_entry

    except Exception as e:
        print(f"✗ Error: {e} (falling back to ChEBI parsed from entry)")
        return chebi_from_entry


def get_chebi_info(chebi_id):
    """Retrieve ChEBI compound information."""
    print(f"\n{'='*70}")
    print("STEP 4: ChEBI Details")
    print(f"{'='*70}")

    if not chebi_id:
        print("⊘ No ChEBI ID available")
        return None

    try:
        c = ChEBI()

        print(f"Retrieving ChEBI entry for {chebi_id}...")

        # Ensure proper format
        if not chebi_id.startswith("CHEBI:"):
            chebi_id = f"CHEBI:{chebi_id}"

        entity = c.getCompleteEntity(chebi_id)

        if entity:
            # ChEBI moved to a REST API in bioservices 1.13; getCompleteEntity now
            # returns a dict-like ChebiEntity whose formula lives on `.formula`
            # (the old SOAP attribute was `Formulae`).
            print(f"\n✓ ChEBI Information:")
            print(f"  ID: {entity.chebiId}")
            print(f"  Name: {entity.chebiAsciiName}")

            if entity.formula:
                print(f"  Formula: {entity.formula}")

            if entity.mass:
                print(f"  Mass: {entity.mass}")

            if entity.charge:
                print(f"  Charge: {entity.charge}")

            return {
                'chebi_id': entity.chebiId,
                'name': entity.chebiAsciiName,
                'formula': entity.formula,
                'mass': entity.mass,
                'smiles': entity.smiles,
            }
        else:
            print("✗ Failed to retrieve ChEBI entry")
            return None

    except Exception as e:
        print(f"✗ Error: {e}")
        return None


def save_results(compound_name, kegg_info, chebi_id, output_file):
    """Save results to file."""
    print(f"\n{'='*70}")
    print(f"Saving results to {output_file}")
    print(f"{'='*70}")

    with open(output_file, 'w') as f:
        f.write("=" * 70 + "\n")
        f.write(f"Compound Cross-Reference Report: {compound_name}\n")
        f.write("=" * 70 + "\n\n")

        # KEGG information
        if kegg_info:
            f.write("KEGG Compound\n")
            f.write("-" * 70 + "\n")
            f.write(f"ID: {kegg_info['kegg_id']}\n")
            f.write(f"Name: {kegg_info['name']}\n")
            f.write(f"Formula: {kegg_info['formula']}\n")
            f.write(f"Exact Mass: {kegg_info['exact_mass']}\n")
            f.write(f"Molecular Weight: {kegg_info['mol_weight']}\n")
            f.write(f"Pathways: {len(kegg_info['pathways'])} found\n")
            f.write("\n")

        # Database IDs
        f.write("Cross-Database Identifiers\n")
        f.write("-" * 70 + "\n")
        if kegg_info:
            f.write(f"KEGG: {kegg_info['kegg_id']}\n")
        if chebi_id:
            f.write(f"ChEBI: {chebi_id}\n")
        f.write("\n")

    print(f"✓ Results saved")


def main():
    """Main workflow."""
    parser = argparse.ArgumentParser(
        description="Search compound across multiple databases",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog="""
Examples:
  python compound_cross_reference.py Geldanamycin
  python compound_cross_reference.py "Adenosine triphosphate"
  python compound_cross_reference.py Aspirin --output aspirin_info.txt
        """
    )
    parser.add_argument("compound", help="Compound name to search")
    parser.add_argument("--output", default=None,
                       help="Output file for results (optional)")

    args = parser.parse_args()

    print("=" * 70)
    print("BIOSERVICES: Compound Cross-Database Search")
    print("=" * 70)

    # Step 1: Search KEGG
    kegg, kegg_id = search_kegg_compound(args.compound)
    if not kegg_id:
        print("\n✗ Failed to find compound. Exiting.")
        sys.exit(1)

    # Step 2: Get KEGG details
    kegg_info = get_kegg_info(kegg, kegg_id)

    # Step 3: Cross-reference KEGG -> ChEBI (via KEGG.conv, with entry fallback)
    chebi_from_entry = kegg_info['chebi_id'] if kegg_info else None
    chebi_id = confirm_chebi_via_conv(kegg, kegg_id, chebi_from_entry)

    # Step 4: Get ChEBI details
    _chebi_info = None
    if chebi_id:
        _chebi_info = get_chebi_info(chebi_id)

    # Summary
    print(f"\n{'='*70}")
    print("SUMMARY")
    print(f"{'='*70}")
    print(f"  Compound: {args.compound}")
    if kegg_info:
        print(f"  KEGG ID: {kegg_info['kegg_id']}")
    if chebi_id:
        print(f"  ChEBI ID: {chebi_id}")
    print("  ChEMBL ID: not mapped here (no bioservices KEGG->ChEMBL route;")
    print("             use the ChEMBL REST API / chembl_webresource_client)")
    print(f"{'='*70}")

    # Save to file if requested
    if args.output:
        save_results(args.compound, kegg_info, chebi_id, args.output)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
KEGG Pathway Network Analysis

This script analyzes all pathways for an organism and extracts:
- Pathway sizes (number of genes)
- Protein-protein interactions
- Interaction type distributions
- Network data in various formats (CSV, SIF)

Usage:
    python pathway_analysis.py ORGANISM OUTPUT_DIR [--limit N]

Examples:
    python pathway_analysis.py hsa ./human_pathways
    python pathway_analysis.py mmu ./mouse_pathways --limit 50

Organism codes:
    hsa = Homo sapiens (human)
    mmu = Mus musculus (mouse)
    dme = Drosophila melanogaster
    sce = Saccharomyces cerevisiae (yeast)
    eco = Escherichia coli
"""

import sys
import os
import argparse
import csv
from collections import Counter
from bioservices import KEGG


def get_all_pathways(kegg, organism):
    """Get all pathway IDs for organism."""
    print(f"\nRetrieving pathways for {organism}...")

    kegg.organism = organism
    pathway_ids = kegg.pathwayIds

    print(f"✓ Found {len(pathway_ids)} pathways")

    return pathway_ids


def analyze_pathway(kegg, pathway_id):
    """Analyze single pathway for size and interactions."""
    try:
        # Parse KGML pathway
        kgml = kegg.parse_kgml_pathway(pathway_id)

        entries = kgml.get('entries', [])
        relations = kgml.get('relations', [])

        # Count relation types
        relation_types = Counter()
        for rel in relations:
            rel_type = rel.get('name', 'unknown')
            relation_types[rel_type] += 1

        # Get pathway name
        try:
            entry = kegg.get(pathway_id)
            pathway_name = "Unknown"
            for line in entry.split("\n"):
                if line.startswith("NAME"):
                    pathway_name = line.replace("NAME", "").strip()
                    break
        except Exception:
            pathway_name = "Unknown"

        result = {
            'pathway_id': pathway_id,
            'pathway_name': pathway_name,
            'num_entries': len(entries),
            'num_relations': len(relations),
            'relation_types': dict(relation_types),
            'entries': entries,
            'relations': relations
        }

        return result

    except Exception as e:
        print(f"  ✗ Error analyzing {pathway_id}: {e}")
        return None


def analyze_all_pathways(kegg, pathway_ids, limit=None):
    """Analyze all pathways."""
    if limit:
        pathway_ids = pathway_ids[:limit]
        print(f"\n⚠ Limiting analysis to first {limit} pathways")

    print(f"\nAnalyzing {len(pathway_ids)} pathways...")

    results = []
    for i, pathway_id in enumerate(pathway_ids, 1):
        print(f"  [{i}/{len(pathway_ids)}] {pathway_id}", end="\r")

        result = analyze_pathway(kegg, pathway_id)
        if result:
            results.append(result)

    print(f"\n✓ Successfully analyzed {len(results)}/{len(pathway_ids)} pathways")

    return results


def save_pathway_summary(results, output_file):
    """Save pathway summary to CSV."""
    print(f"\nSaving pathway summary to {output_file}...")

    with open(output_file, 'w', newline='') as f:
        writer = csv.writer(f)

        # Header
        writer.writerow([
            'Pathway_ID',
            'Pathway_Name',
            'Num_Genes',
            'Num_Interactions',
            'Activation',
            'Inhibition',
            'Phosphorylation',
            'Binding',
            'Other'
        ])

        # Data
        for result in results:
            rel_types = result['relation_types']

            writer.writerow([
                result['pathway_id'],
                result['pathway_name'],
                result['num_entries'],
                result['num_relations'],
                rel_types.get('activation', 0),
                rel_types.get('inhibition', 0),
                rel_types.get('phosphorylation', 0),
                rel_types.get('binding/association', 0),
                sum(v for k, v in rel_types.items()
                    if k not in ['activation', 'inhibition', 'phosphorylation', 'binding/association'])
            ])

    print(f"✓ Summary saved")


def save_interactions_sif(results, output_file):
    """Save all interactions in SIF format."""
    print(f"\nSaving interactions to {output_file}...")

    with open(output_file, 'w') as f:
        for result in results:
            _pathway_id = result['pathway_id']
            for rel in result['relations']:
                entry1 = rel.get('entry1', '')
                entry2 = rel.get('entry2', '')
                interaction_type = rel.get('name', 'interaction')

                # Write SIF format: source\tinteraction\ttarget
                f.write(f"{entry1}\t{interaction_type}\t{entry2}\n")

    print(f"✓ Interactions saved")


def save_detailed_pathway_info(results, output_dir):
    """Save detailed information for each pathway."""
    print(f"\nSaving detailed pathway files to {output_dir}/pathways/...")

    pathway_dir = os.path.join(output_dir, "pathways")
    os.makedirs(pathway_dir, exist_ok=True)

    for result in results:
        pathway_id = result['pathway_id'].replace(":", "_")
        filename = os.path.join(pathway_dir, f"{pathway_id}_interactions.csv")

        with open(filename, 'w', newline='') as f:
            writer = csv.writer(f)
            writer.writerow(['Source', 'Target', 'Interaction_Type', 'Link_Type'])

            for rel in result['relations']:
                writer.writerow([
                    rel.get('entry1', ''),
                    rel.get('entry2', ''),
                    rel.get('name', 'unknown'),
                    rel.get('link', 'unknown')
                ])

    print(f"✓ Detailed files saved for {len(results)} pathways")


def print_statistics(results):
    """Print analysis statistics."""
    print(f"\n{'='*70}")
    print("PATHWAY ANALYSIS STATISTICS")
    print(f"{'='*70}")

    # Total stats
    total_pathways = len(results)
    total_interactions = sum(r['num_relations'] for r in results)
    total_genes = sum(r['num_entries'] for r in results)

    print(f"\nOverall:")
    print(f"  Total pathways: {total_pathways}")
    print(f"  Total genes/proteins: {total_genes}")
    print(f"  Total interactions: {total_interactions}")

    # Largest pathways
    print(f"\nLargest pathways (by gene count):")
    sorted_by_size = sorted(results, key=lambda x: x['num_entries'], reverse=True)
    for i, result in enumerate(sorted_by_size[:10], 1):
        print(f"  {i}. {result['pathway_id']}: {result['num_entries']} genes")
        print(f"     {result['pathway_name']}")

    # Most connected pathways
    print(f"\nMost connected pathways (by interactions):")
    sorted_by_connections = sorted(results, key=lambda x: x['num_relations'], reverse=True)
    for i, result in enumerate(sorted_by_connections[:10], 1):
        print(f"  {i}. {result['pathway_id']}: {result['num_relations']} interactions")
        print(f"     {result['pathway_name']}")

    # Interaction type distribution
    print(f"\nInteraction type distribution:")
    all_types = Counter()
    for result in results:
        for rel_type, count in result['relation_types'].items():
            all_types[rel_type] += count

    for rel_type, count in all_types.most_common():
        percentage = (count / total_interactions) * 100 if total_interactions > 0 else 0
        print(f"  {rel_type}: {count} ({percentage:.1f}%)")


def main():
    """Main analysis workflow."""
    parser = argparse.ArgumentParser(
        description="Analyze KEGG pathways for an organism",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog="""
Examples:
  python pathway_analysis.py hsa ./human_pathways
  python pathway_analysis.py mmu ./mouse_pathways --limit 50

Organism codes:
  hsa = Homo sapiens (human)
  mmu = Mus musculus (mouse)
  dme = Drosophila melanogaster
  sce = Saccharomyces cerevisiae (yeast)
  eco = Escherichia coli
        """
    )
    parser.add_argument("organism", help="KEGG organism code (e.g., hsa, mmu)")
    parser.add_argument("output_dir", help="Output directory for results")
    parser.add_argument("--limit", type=int, default=None,
                       help="Limit analysis to first N pathways")

    args = parser.parse_args()

    print("=" * 70)
    print("BIOSERVICES: KEGG Pathway Network Analysis")
    print("=" * 70)

    # Create output directory
    os.makedirs(args.output_dir, exist_ok=True)

    # Initialize KEGG
    kegg = KEGG()

    # Get all pathways
    pathway_ids = get_all_pathways(kegg, args.organism)

    if not pathway_ids:
        print(f"\n✗ No pathways found for {args.organism}")
        sys.exit(1)

    # Analyze pathways
    results = analyze_all_pathways(kegg, pathway_ids, args.limit)

    if not results:
        print("\n✗ No pathways successfully analyzed")
        sys.exit(1)

    # Print statistics
    print_statistics(results)

    # Save results
    summary_file = os.path.join(args.output_dir, "pathway_summary.csv")
    save_pathway_summary(results, summary_file)

    sif_file = os.path.join(args.output_dir, "all_interactions.sif")
    save_interactions_sif(results, sif_file)

    save_detailed_pathway_info(results, args.output_dir)

    # Final summary
    print(f"\n{'='*70}")
    print("OUTPUT FILES")
    print(f"{'='*70}")
    print(f"  Summary: {summary_file}")
    print(f"  Interactions: {sif_file}")
    print(f"  Detailed: {args.output_dir}/pathways/")
    print(f"{'='*70}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Complete Protein Analysis Workflow

This script performs a comprehensive protein analysis pipeline:
1. UniProt search and identifier retrieval
2. FASTA sequence retrieval
3. BLAST similarity search
4. KEGG pathway discovery
5. STRING interaction mapping
6. GO annotation retrieval

Usage:
    python protein_analysis_workflow.py PROTEIN_NAME EMAIL [--skip-blast]

Examples:
    python protein_analysis_workflow.py ZAP70_HUMAN user@example.com
    python protein_analysis_workflow.py P43403 user@example.com --skip-blast

Note: BLAST searches can take several minutes. Use --skip-blast to skip this step.
"""

import sys
import time
import argparse
from bioservices import UniProt, KEGG, NCBIblast, STRING, QuickGO


def search_protein(query):
    """Search UniProt for protein and retrieve basic information."""
    print(f"\n{'='*70}")
    print("STEP 1: UniProt Search")
    print(f"{'='*70}")

    u = UniProt(verbose=False)

    print(f"Searching for: {query}")

    # Try direct retrieval first (if query looks like accession)
    if len(query) == 6 and query[0] in "OPQ":
        try:
            entry = u.retrieve(query, frmt="txt")
            if entry:
                uniprot_id = query
                print(f"✓ Found UniProt entry: {uniprot_id}")
                return u, uniprot_id
        except Exception:
            pass

    # Otherwise search
    results = u.search(
        query,
        frmt="tsv",
        columns="accession,gene_names,organism_name,length,protein_name",
        limit=5,
    )

    if not results:
        print("✗ No results found")
        return u, None

    lines = results.strip().split("\n")
    if len(lines) < 2:
        print("✗ No entries found")
        return u, None

    # Display results
    print(f"\n✓ Found {len(lines)-1} result(s):")
    for i, line in enumerate(lines[1:], 1):
        fields = line.split("\t")
        print(f"  {i}. {fields[0]} - {fields[1]} ({fields[2]})")

    # Use first result
    first_entry = lines[1].split("\t")
    uniprot_id = first_entry[0]
    gene_names = first_entry[1] if len(first_entry) > 1 else "N/A"
    organism = first_entry[2] if len(first_entry) > 2 else "N/A"
    length = first_entry[3] if len(first_entry) > 3 else "N/A"
    protein_name = first_entry[4] if len(first_entry) > 4 else "N/A"

    print(f"\nUsing first result:")
    print(f"  UniProt ID: {uniprot_id}")
    print(f"  Gene names: {gene_names}")
    print(f"  Organism: {organism}")
    print(f"  Length: {length} aa")
    print(f"  Protein: {protein_name}")

    return u, uniprot_id


def retrieve_sequence(uniprot, uniprot_id):
    """Retrieve FASTA sequence for protein."""
    print(f"\n{'='*70}")
    print("STEP 2: FASTA Sequence Retrieval")
    print(f"{'='*70}")

    try:
        sequence = uniprot.retrieve(uniprot_id, frmt="fasta")

        if sequence:
            # Extract sequence only (remove header)
            lines = sequence.strip().split("\n")
            header = lines[0]
            seq_only = "".join(lines[1:])

            print(f"✓ Retrieved sequence:")
            print(f"  Header: {header}")
            print(f"  Length: {len(seq_only)} residues")
            print(f"  First 60 residues: {seq_only[:60]}...")

            return seq_only
        else:
            print("✗ Failed to retrieve sequence")
            return None

    except Exception as e:
        print(f"✗ Error: {e}")
        return None


def run_blast(sequence, email, skip=False):
    """Run BLAST similarity search."""
    print(f"\n{'='*70}")
    print("STEP 3: BLAST Similarity Search")
    print(f"{'='*70}")

    if skip:
        print("⊘ Skipped (--skip-blast flag)")
        return None

    if not email or "@" not in email:
        print("⊘ Skipped (valid email required for BLAST)")
        return None

    try:
        print(f"Submitting BLASTP job...")
        print(f"  Database: uniprotkb")
        print(f"  Sequence length: {len(sequence)} aa")

        s = NCBIblast(verbose=False)

        jobid = s.run(
            program="blastp",
            sequence=sequence,
            stype="protein",
            database="uniprotkb",
            email=email
        )

        print(f"✓ Job submitted: {jobid}")
        print(f"  Waiting for completion...")

        # Poll for completion
        max_wait = 300  # 5 minutes
        start_time = time.time()

        while time.time() - start_time < max_wait:
            status = s.get_status(jobid)
            elapsed = int(time.time() - start_time)
            print(f"  Status: {status} (elapsed: {elapsed}s)", end="\r")

            if status == "FINISHED":
                print(f"\n✓ BLAST completed in {elapsed}s")

                # Retrieve results
                results = s.get_result(jobid, "out")

                # Parse and display summary
                lines = results.split("\n")
                print(f"\n  Results preview:")
                for line in lines[:20]:
                    if line.strip():
                        print(f"    {line}")

                return results

            elif status in ("ERROR", "FAILURE", "NOT_FOUND"):
                print(f"\n✗ BLAST job did not complete: {status}")
                return None

            time.sleep(5)

        print(f"\n✗ Timeout after {max_wait}s")
        return None

    except Exception as e:
        print(f"✗ Error: {e}")
        return None


def discover_pathways(uniprot, kegg, uniprot_id):
    """Discover KEGG pathways for protein."""
    print(f"\n{'='*70}")
    print("STEP 4: KEGG Pathway Discovery")
    print(f"{'='*70}")

    try:
        # Map UniProt → KEGG. mapping() returns UniProt's job payload:
        # {"results": [{"from": ..., "to": ...}], "failedIds": [...]}
        print(f"Mapping {uniprot_id} to KEGG...")
        job = uniprot.mapping(fr="UniProtKB_AC-ID", to="KEGG", query=uniprot_id)

        kegg_ids = [
            row["to"] for row in (job or {}).get("results", []) if row["from"] == uniprot_id
        ]
        if not kegg_ids:
            print("✗ No KEGG mapping found")
            return []
        print(f"✓ KEGG ID(s): {kegg_ids}")

        # Get pathways for first KEGG ID
        kegg_id = kegg_ids[0]
        organism, gene_id = kegg_id.split(":")

        print(f"\nSearching pathways for {kegg_id}...")
        pathways = kegg.get_pathway_by_gene(gene_id, organism)

        if not pathways:
            print("✗ No pathways found")
            return []

        print(f"✓ Found {len(pathways)} pathway(s):\n")

        # get_pathway_by_gene returns the parsed PATHWAY block, which current
        # bioservices gives as {pathway_id: name}; tolerate a plain sequence too.
        if isinstance(pathways, dict):
            pathway_pairs = list(pathways.items())
        else:
            pathway_pairs = [(pid, None) for pid in pathways]

        pathway_info = []
        for pathway_id, pathway_name in pathway_pairs:
            if not pathway_name:
                try:
                    entry = kegg.get(pathway_id)
                    pathway_name = "Unknown"
                    for line in entry.split("\n"):
                        if line.startswith("NAME"):
                            pathway_name = line.replace("NAME", "").strip()
                            break
                except Exception:
                    pathway_name = "[Error retrieving name]"

            pathway_info.append((pathway_id, pathway_name))
            print(f"  • {pathway_id}: {pathway_name}")

        return pathway_info

    except Exception as e:
        print(f"✗ Error: {e}")
        return []


def find_interactions(protein_query):
    """Find protein-protein interactions via STRING."""
    print(f"\n{'='*70}")
    print("STEP 5: Protein-Protein Interactions")
    print(f"{'='*70}")

    try:
        s = STRING()

        # PSICQUIC and BioGRID were removed from bioservices in 1.14; STRING answers
        # the same question. species is an NCBI taxid, required_score is 0-1000.
        print("Querying STRING for human interaction partners...")
        print(f"  Query: {protein_query} (taxid 9606, score >= 400)")

        rows = s.get_interaction_partners(
            protein_query, species=9606, required_score=400, limit=25
        )

        if not rows:
            print("✗ No interactions found in STRING")
            return []

        print(f"✓ Found {len(rows)} interaction(s):\n")

        interactions = []
        for i, row in enumerate(rows[:10], 1):
            protein_a = row.get("preferredName_A", "?")
            protein_b = row.get("preferredName_B", "?")
            score = row.get("score")
            interactions.append((protein_a, protein_b, score))
            print(f"  {i}. {protein_a} ↔ {protein_b} (score {score})")

        if len(rows) > 10:
            print(f"  ... and {len(rows)-10} more")

        return interactions

    except Exception as e:
        print(f"✗ Error: {e}")
        return []


def get_go_annotations(uniprot_id):
    """Retrieve GO annotations."""
    print(f"\n{'='*70}")
    print("STEP 6: Gene Ontology Annotations")
    print(f"{'='*70}")

    try:
        g = QuickGO()

        print(f"Retrieving GO annotations for {uniprot_id}...")
        # QuickGO REST parameters: geneProductId is prefixed and limit maxes out at 100.
        annotations = g.Annotation(
            geneProductId=f"UniProtKB:{uniprot_id}",
            includeFields="goName",
            limit=100,
        )

        rows = (annotations or {}).get("results") if isinstance(annotations, dict) else None
        if not rows:
            print("✗ No GO annotations found")
            return {}

        print(f"✓ Found {annotations.get('numberOfHits', len(rows))} annotation(s)\n")

        # Group by aspect; goAspect is spelled out (biological_process, ...)
        aspect_keys = {"biological_process": "P", "molecular_function": "F", "cellular_component": "C"}
        aspects = {"P": [], "F": [], "C": []}
        for row in rows:
            key = aspect_keys.get(row.get("goAspect"))
            if key:
                aspects[key].append((row["goId"], row.get("goName", "")))

        # Display summary
        print(f"  Biological Process (P): {len(aspects['P'])} terms")
        for go_id, go_term in aspects['P'][:5]:
            print(f"    • {go_id}: {go_term}")
        if len(aspects['P']) > 5:
            print(f"    ... and {len(aspects['P'])-5} more")

        print(f"\n  Molecular Function (F): {len(aspects['F'])} terms")
        for go_id, go_term in aspects['F'][:5]:
            print(f"    • {go_id}: {go_term}")
        if len(aspects['F']) > 5:
            print(f"    ... and {len(aspects['F'])-5} more")

        print(f"\n  Cellular Component (C): {len(aspects['C'])} terms")
        for go_id, go_term in aspects['C'][:5]:
            print(f"    • {go_id}: {go_term}")
        if len(aspects['C']) > 5:
            print(f"    ... and {len(aspects['C'])-5} more")

        return aspects

    except Exception as e:
        print(f"✗ Error: {e}")
        return {}


def main():
    """Main workflow."""
    parser = argparse.ArgumentParser(
        description="Complete protein analysis workflow using BioServices",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog="""
Examples:
  python protein_analysis_workflow.py ZAP70_HUMAN user@example.com
  python protein_analysis_workflow.py P43403 user@example.com --skip-blast
        """
    )
    parser.add_argument("protein", help="Protein name or UniProt ID")
    parser.add_argument("email", help="Email address (required for BLAST)")
    parser.add_argument("--skip-blast", action="store_true",
                       help="Skip BLAST search (faster)")

    args = parser.parse_args()

    print("=" * 70)
    print("BIOSERVICES: Complete Protein Analysis Workflow")
    print("=" * 70)

    # Step 1: Search protein
    uniprot, uniprot_id = search_protein(args.protein)
    if not uniprot_id:
        print("\n✗ Failed to find protein. Exiting.")
        sys.exit(1)

    # Step 2: Retrieve sequence
    sequence = retrieve_sequence(uniprot, uniprot_id)
    if not sequence:
        print("\n⚠ Warning: Could not retrieve sequence")

    # Step 3: BLAST search
    if sequence:
        _blast_results = run_blast(sequence, args.email, args.skip_blast)
    # Step 4: Pathway discovery
    kegg = KEGG()
    pathways = discover_pathways(uniprot, kegg, uniprot_id)

    # Step 5: Interaction mapping
    interactions = find_interactions(args.protein)

    # Step 6: GO annotations
    go_terms = get_go_annotations(uniprot_id)

    # Summary
    print(f"\n{'='*70}")
    print("WORKFLOW SUMMARY")
    print(f"{'='*70}")
    print(f"  Protein: {args.protein}")
    print(f"  UniProt ID: {uniprot_id}")
    print(f"  Sequence: {'✓' if sequence else '✗'}")
    print(f"  BLAST: {'✓' if not args.skip_blast and sequence else '⊘'}")
    print(f"  Pathways: {len(pathways)} found")
    print(f"  Interactions: {len(interactions)} found")
    print(f"  GO annotations: {sum(len(v) for v in go_terms.values())} found")
    print(f"{'='*70}")


if __name__ == "__main__":
    main()