Skills › Scientific Agent Skills › bioservices › protein_analysis_workflow.py
protein_analysis_workflow.py
377 lines · skills/bioservices/scripts/protein_analysis_workflow.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Complete Protein Analysis Workflow
5This script performs a comprehensive protein analysis pipeline:
61. UniProt search and identifier retrieval
72. FASTA sequence retrieval
83. BLAST similarity search
94. KEGG pathway discovery
105. STRING association lookup
116. GO annotation retrieval
13Usage:
14 export NCBI_EMAIL=you@lab.org
15 python protein_analysis_workflow.py PROTEIN_NAME [EMAIL] [--skip-blast]
17Examples:
18 python protein_analysis_workflow.py ZAP70_HUMAN
19 python protein_analysis_workflow.py P43403 user@example.com --skip-blast
21Note: BLAST searches can take several minutes. Use --skip-blast to skip this step.
22Email is read from NCBI_EMAIL when the optional EMAIL argument is omitted.
23"""
25import os
26import re
27import sys
28import time
29import argparse
30from bioservices import UniProt, KEGG, NCBIblast, QuickGO, STRING
31from batch_id_converter import mapping_to_lists
33_EMAIL_RE = re.compile(r"^[^@\s]+@[^@\s]+\.[^@\s]+$")
36def resolve_ncbi_email(cli_email=None):
37 """Return a validated EMBL-EBI BLAST contact email from CLI or NCBI_EMAIL."""
38 email = (cli_email or os.environ.get("NCBI_EMAIL", "")).strip()
39 if email and _EMAIL_RE.match(email):
40 return email
41 return None
44def search_protein(query):
45 """Search UniProt for protein and retrieve basic information."""
46 print(f"\n{'='*70}")
47 print("STEP 1: UniProt Search")
48 print(f"{'='*70}")
50 u = UniProt(verbose=False)
52 print(f"Searching for: {query}")
54 # Try direct retrieval first (if query looks like accession)
55 if re.fullmatch(r"[A-Z0-9]{6}|[A-Z0-9]{10}", query):
56 try:
57 entry = u.retrieve(query, frmt="json")
58 if isinstance(entry, dict) and entry.get("primaryAccession"):
59 uniprot_id = entry["primaryAccession"]
60 print(f"[OK] Found UniProt entry: {uniprot_id}")
61 return u, uniprot_id
62 except:
63 pass
65 # Otherwise search
66 results = u.search(query, frmt="tsv", columns="accession,gene_names,organism_name,length,protein_name", limit=5, size=5)
68 if not isinstance(results, str) or not results.strip():
69 print("[FAIL] Search failed or returned no results")
70 return u, None
72 lines = results.strip().split("\n")
73 if len(lines) < 2:
74 print("[FAIL] No entries found")
75 return u, None
77 # Display results
78 print(f"\n[OK] Found {len(lines)-1} result(s):")
79 for i, line in enumerate(lines[1:], 1):
80 fields = line.split("\t")
81 print(f" {i}. {fields[0]} - {fields[1]} ({fields[2]})")
83 # Use first result
84 first_entry = lines[1].split("\t")
85 uniprot_id = first_entry[0]
86 gene_names = first_entry[1] if len(first_entry) > 1 else "N/A"
87 organism = first_entry[2] if len(first_entry) > 2 else "N/A"
88 length = first_entry[3] if len(first_entry) > 3 else "N/A"
89 protein_name = first_entry[4] if len(first_entry) > 4 else "N/A"
91 print(f"\nUsing first result:")
92 print(f" UniProt ID: {uniprot_id}")
93 print(f" Gene names: {gene_names}")
94 print(f" Organism: {organism}")
95 print(f" Length: {length} aa")
96 print(f" Protein: {protein_name}")
98 return u, uniprot_id
101def retrieve_sequence(uniprot, uniprot_id):
102 """Retrieve FASTA sequence for protein."""
103 print(f"\n{'='*70}")
104 print("STEP 2: FASTA Sequence Retrieval")
105 print(f"{'='*70}")
107 try:
108 sequence = uniprot.retrieve(uniprot_id, frmt="fasta")
110 if isinstance(sequence, str) and sequence.startswith(">"):
111 # Extract sequence only (remove header)
112 lines = sequence.strip().split("\n")
113 header = lines[0]
114 seq_only = "".join(lines[1:])
116 print(f"[OK] Retrieved sequence:")
117 print(f" Header: {header}")
118 print(f" Length: {len(seq_only)} residues")
119 print(f" First 60 residues: {seq_only[:60]}...")
121 return seq_only
122 else:
123 print("[FAIL] Failed to retrieve sequence")
124 return None
126 except Exception as e:
127 print(f"[FAIL] Error: {e}")
128 return None
131def run_blast(sequence, email, skip=False):
132 """Run BLAST similarity search."""
133 print(f"\n{'='*70}")
134 print("STEP 3: BLAST Similarity Search")
135 print(f"{'='*70}")
137 if skip:
138 print("[SKIP] Skipped (--skip-blast flag)")
139 return None
141 if not email:
142 print("[SKIP] Skipped (set NCBI_EMAIL or pass email for BLAST)")
143 return None
145 try:
146 print(f"Submitting BLASTP job...")
147 print(f" Database: uniprotkb")
148 print(f" Sequence length: {len(sequence)} aa")
150 s = NCBIblast(verbose=False)
151 s.services.url = "https://www.ebi.ac.uk/Tools/services/rest/ncbiblast"
153 jobid = s.run(
154 program="blastp",
155 sequence=sequence,
156 stype="protein",
157 database="uniprotkb",
158 email=email
159 )
161 if not isinstance(jobid, str) or not jobid.startswith("ncbiblast-"):
162 raise ValueError(f"Invalid BLAST job response: {jobid!r}")
163 print(f"[OK] Job submitted: {jobid}")
164 print(f" Waiting for completion...")
166 # Poll for completion
167 max_wait = 300 # 5 minutes
168 start_time = time.time()
170 while time.time() - start_time < max_wait:
171 status = s.get_status(jobid)
172 elapsed = int(time.time() - start_time)
173 print(f" Status: {status} (elapsed: {elapsed}s)", end="\r")
175 if status == "FINISHED":
176 print(f"\n[OK] BLAST completed in {elapsed}s")
178 # Retrieve results
179 results = s.get_result(jobid, "out")
181 # Parse and display summary
182 lines = results.split("\n")
183 print(f"\n Results preview:")
184 for line in lines[:20]:
185 if line.strip():
186 print(f" {line}")
188 return results
190 elif status in {"ERROR", "FAILURE", "NOT_FOUND"}:
191 print(f"\n[FAIL] BLAST job failed")
192 return None
194 time.sleep(5)
196 print(f"\n[FAIL] Timeout after {max_wait}s")
197 return None
199 except Exception as e:
200 print(f"[FAIL] Error: {e}")
201 return None
204def discover_pathways(uniprot, kegg, uniprot_id):
205 """Discover KEGG pathways for protein."""
206 print(f"\n{'='*70}")
207 print("STEP 4: KEGG Pathway Discovery")
208 print(f"{'='*70}")
210 try:
211 # Map UniProt → KEGG
212 print(f"Mapping {uniprot_id} to KEGG...")
213 kegg_mapping = mapping_to_lists(uniprot.mapping(
214 fr="UniProtKB_AC-ID", to="KEGG", query=uniprot_id))
215 kegg_ids = kegg_mapping.get(uniprot_id) or []
216 if not kegg_ids:
217 print("[SKIP] No KEGG mapping returned")
218 return []
219 pathway_info = {}
220 for kegg_id in kegg_ids:
221 organism, gene_id = kegg_id.split(":", 1)
222 found = kegg.get_pathway_by_gene(gene_id, organism)
223 if found is None:
224 continue
225 if not isinstance(found, dict):
226 raise ValueError("Unexpected KEGG PATHWAY response")
227 pathway_info.update(found)
228 for pathway_id, name in pathway_info.items():
229 print(f" {pathway_id}: {name}")
230 return list(pathway_info.items())
232 except Exception as e:
233 print(f"[FAIL] Error: {e}")
234 return []
237def find_interactions(protein_query, taxon_id=None):
238 """Return up to ten STRING associations, retaining scores and stable IDs."""
239 if taxon_id is None:
240 print("[SKIP] STRING lookup needs the protein's verified taxonomy ID")
241 return []
242 try:
243 results = STRING(verbose=False).get_interaction_partners(
244 protein_query, species=taxon_id, required_score=700, limit=10,
245 caller_identity="scientific-agent-skills-bioservices")
246 if not isinstance(results, list):
247 raise ValueError("Unexpected STRING response")
248 print(f"[OK] {len(results)} STRING associations (not necessarily physical binding)")
249 return results
250 except Exception as error:
251 print(f"[FAIL] STRING lookup: {error}")
252 return []
255def get_go_annotations(uniprot_id, max_pages=1000):
256 """Fetch all QuickGO pages and group distinct positive terms by aspect.
258 NOT-qualified annotations are not positive assertions. Evidence and complete
259 raw records should be retained separately for downstream enrichment.
260 """
261 if type(max_pages) is not int or max_pages < 1:
262 raise ValueError("max_pages must be a positive integer")
263 g = QuickGO(verbose=False)
264 aspects = {"P": [], "F": [], "C": []}
265 aspect_codes = {"biological_process": "P", "molecular_function": "F",
266 "cellular_component": "C"}
267 expected_total = None
268 for page in range(1, max_pages + 1):
269 payload = g.Annotation(geneProductId=f"UniProtKB:{uniprot_id}",
270 includeFields="goName", limit=100, page=page)
271 if not isinstance(payload, dict) or not isinstance(payload.get("results"), list):
272 raise ValueError("QuickGO lookup failed; annotations are incomplete")
273 page_info = payload.get("pageInfo")
274 if not isinstance(page_info, dict):
275 raise ValueError("QuickGO response lacks pagination metadata")
276 total_pages = page_info.get("total")
277 if type(total_pages) is not int or not 0 <= total_pages <= max_pages:
278 raise ValueError("QuickGO page total is invalid or exceeds max_pages; annotations are incomplete")
279 current = page_info.get("current", page)
280 if type(current) is not int or current != page:
281 raise ValueError("QuickGO returned the wrong page; annotations are incomplete")
282 if expected_total is not None and total_pages != expected_total:
283 raise ValueError("QuickGO page total changed; rerun to obtain complete annotations")
284 expected_total = total_pages
285 if total_pages == 0 and payload["results"]:
286 raise ValueError("QuickGO zero-page response contains results")
287 for record in payload["results"]:
288 if "NOT" in str(record.get("qualifier", "")).split("|"):
289 continue
290 aspect = aspect_codes.get(record.get("goAspect"))
291 term = (record.get("goId"), record.get("goName", ""))
292 if aspect and term[0] and term not in aspects[aspect]:
293 aspects[aspect].append(term)
294 if page >= total_pages:
295 break
296 print(f"[OK] {sum(map(len, aspects.values()))} distinct positive GO terms")
297 return aspects
300def main():
301 """Main workflow."""
302 parser = argparse.ArgumentParser(
303 description="Complete protein analysis workflow using BioServices",
304 formatter_class=argparse.RawDescriptionHelpFormatter,
305 epilog="""
306Examples:
307 export NCBI_EMAIL=you@lab.org
308 python protein_analysis_workflow.py ZAP70_HUMAN
309 python protein_analysis_workflow.py P43403 user@example.com --skip-blast
310 """
311 )
312 parser.add_argument("protein", help="Protein name or UniProt ID")
313 parser.add_argument(
314 "email",
315 nargs="?",
316 default=None,
317 help="EMBL-EBI BLAST contact email (optional if NCBI_EMAIL is set)",
318 )
319 parser.add_argument("--skip-blast", action="store_true",
320 help="Skip BLAST search (faster)")
322 args = parser.parse_args()
324 print("=" * 70)
325 print("BIOSERVICES: Complete Protein Analysis Workflow")
326 print("=" * 70)
328 # Step 1: Search protein
329 uniprot, uniprot_id = search_protein(args.protein)
330 if not uniprot_id:
331 print("\n[FAIL] Failed to find protein. Exiting.")
332 sys.exit(1)
334 # Step 2: Retrieve sequence
335 sequence = retrieve_sequence(uniprot, uniprot_id)
336 if not sequence:
337 print("\n[WARN] Warning: Could not retrieve sequence")
339 # Step 3: BLAST search
340 ncbi_email = resolve_ncbi_email(args.email)
341 blast_results = None
342 if sequence:
343 blast_results = run_blast(sequence, ncbi_email, args.skip_blast)
345 # Step 4: Pathway discovery
346 kegg = KEGG()
347 kegg.services.url = "https://rest.kegg.jp"
348 pathways = discover_pathways(uniprot, kegg, uniprot_id)
350 # Step 5: Interaction mapping
351 entry = uniprot.retrieve(uniprot_id, frmt="json")
352 taxon_id = entry.get("organism", {}).get("taxonId") if isinstance(entry, dict) else None
353 interactions = find_interactions(uniprot_id, taxon_id)
355 # Step 6: GO annotations
356 try:
357 go_terms = get_go_annotations(uniprot_id)
358 except Exception as error:
359 print(f"[FAIL] {error}")
360 go_terms = None
362 # Summary
363 print(f"\n{'='*70}")
364 print("WORKFLOW SUMMARY")
365 print(f"{'='*70}")
366 print(f" Protein: {args.protein}")
367 print(f" UniProt ID: {uniprot_id}")
368 print(f" Sequence: {'[OK]' if sequence else '[FAIL]'}")
369 print(f" BLAST: {'[OK]' if blast_results else '[SKIP/FAIL]'}")
370 print(f" Pathways: {len(pathways)} found")
371 print(f" Interactions: {len(interactions)} found")
372 print(f" GO terms: {sum(map(len, go_terms.values())) if go_terms is not None else 'unavailable'}")
373 print(f"{'='*70}")
376if __name__ == "__main__":
377 main()