Skills › Scientific Agent Skills › bioservices › pathway_analysis.py
pathway_analysis.py
313 lines · skills/bioservices/scripts/pathway_analysis.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3KEGG Pathway Network Analysis
5This script analyzes all pathways for an organism and extracts:
6- Pathway sizes (number of KGML entries)
7- KGML relation subtype records
8- Interaction type distributions
9- Network data in various formats (CSV, SIF)
11Usage:
12 python pathway_analysis.py ORGANISM OUTPUT_DIR [--limit N]
14Examples:
15 python pathway_analysis.py hsa ./human_pathways
16 python pathway_analysis.py mmu ./mouse_pathways --limit 50
18Organism codes:
19 hsa = Homo sapiens (human)
20 mmu = Mus musculus (mouse)
21 dme = Drosophila melanogaster
22 sce = Saccharomyces cerevisiae (yeast)
23 eco = Escherichia coli
24"""
26import sys
27import os
28import argparse
29import csv
30from urllib.parse import quote
31from collections import Counter
32from bioservices import KEGG
35def get_all_pathways(kegg, organism):
36 """Get all pathway IDs for organism."""
37 print(f"\nRetrieving pathways for {organism}...")
39 raw = kegg.services.http_get("list/pathway/" + quote(organism, safe=""), frmt="txt")
40 if not isinstance(raw, str) or (raw.strip() and "\t" not in raw):
41 raise ValueError("KEGG pathway listing failed or returned an unexpected response")
42 pathway_ids = [line.split("\t", 1)[0] for line in raw.splitlines() if line]
44 print(f"[OK] Found {len(pathway_ids)} pathways")
46 return pathway_ids
49def analyze_pathway(kegg, pathway_id):
50 """Analyze single pathway for size and interactions."""
51 try:
52 # Parse KGML pathway
53 kgml = kegg.parse_kgml_pathway(pathway_id)
55 entries = kgml.get('entries', [])
56 relations = kgml.get('relations', [])
58 # Count relation types
59 relation_types = Counter()
60 for rel in relations:
61 rel_type = rel.get('name') or 'unknown'
62 relation_types[rel_type] += 1
64 # Get pathway name
65 try:
66 entry = kegg.get(pathway_id)
67 pathway_name = "Unknown"
68 for line in entry.split("\n"):
69 if line.startswith("NAME"):
70 pathway_name = line.replace("NAME", "").strip()
71 break
72 except:
73 pathway_name = "Unknown"
75 result = {
76 'pathway_id': pathway_id,
77 'pathway_name': pathway_name,
78 'num_entries': len(entries),
79 'num_relations': len(relations),
80 'relation_types': dict(relation_types),
81 'entries': entries,
82 'relations': relations
83 }
85 return result
87 except Exception as e:
88 print(f" [FAIL] Error analyzing {pathway_id}: {e}")
89 return None
92def analyze_all_pathways(kegg, pathway_ids, limit=None):
93 """Analyze all pathways."""
94 if limit:
95 pathway_ids = pathway_ids[:limit]
96 print(f"\n[WARN] Limiting analysis to first {limit} pathways")
98 print(f"\nAnalyzing {len(pathway_ids)} pathways...")
100 results = []
101 for i, pathway_id in enumerate(pathway_ids, 1):
102 print(f" [{i}/{len(pathway_ids)}] {pathway_id}", end="\r")
104 result = analyze_pathway(kegg, pathway_id)
105 if result:
106 results.append(result)
108 print(f"\n[OK] Successfully analyzed {len(results)}/{len(pathway_ids)} pathways")
110 return results
113def save_pathway_summary(results, output_file):
114 """Save pathway summary to CSV."""
115 print(f"\nSaving pathway summary to {output_file}...")
117 with open(output_file, 'w', newline='') as f:
118 writer = csv.writer(f)
120 # Header
121 writer.writerow([
122 'Pathway_ID',
123 'Pathway_Name',
124 'Num_Entries',
125 'Num_Relation_Records',
126 'Activation',
127 'Inhibition',
128 'Phosphorylation',
129 'Binding',
130 'Other'
131 ])
133 # Data
134 for result in results:
135 rel_types = result['relation_types']
137 writer.writerow([
138 result['pathway_id'],
139 result['pathway_name'],
140 result['num_entries'],
141 result['num_relations'],
142 rel_types.get('activation', 0),
143 rel_types.get('inhibition', 0),
144 rel_types.get('phosphorylation', 0),
145 rel_types.get('binding/association', 0),
146 sum(v for k, v in rel_types.items()
147 if k not in ['activation', 'inhibition', 'phosphorylation', 'binding/association'])
148 ])
150 print(f"[OK] Summary saved")
153def save_interactions_sif(results, output_file):
154 """Save all interactions in SIF format."""
155 print(f"\nSaving interactions to {output_file}...")
157 with open(output_file, 'w') as f:
158 for result in results:
159 pathway_id = result['pathway_id']
161 for rel in result['relations']:
162 entry1 = rel.get('entry1', '')
163 entry2 = rel.get('entry2', '')
164 interaction_type = rel.get('name') or 'interaction'
166 # Write SIF format: source\tinteraction\ttarget
167 f.write(f"{pathway_id}#{entry1}\t{interaction_type}\t{pathway_id}#{entry2}\n")
169 print(f"[OK] Interactions saved")
172def save_detailed_pathway_info(results, output_dir):
173 """Save detailed information for each pathway."""
174 print(f"\nSaving detailed pathway files to {output_dir}/pathways/...")
176 pathway_dir = os.path.join(output_dir, "pathways")
177 os.makedirs(pathway_dir, exist_ok=True)
179 for result in results:
180 pathway_id = result['pathway_id'].replace(":", "_")
181 filename = os.path.join(pathway_dir, f"{pathway_id}_interactions.csv")
183 with open(filename, 'w', newline='') as f:
184 writer = csv.writer(f)
185 writer.writerow(['Source', 'Target', 'Interaction_Type', 'Link_Type'])
187 for rel in result['relations']:
188 writer.writerow([
189 rel.get('entry1', ''),
190 rel.get('entry2', ''),
191 rel.get('name') or 'unknown',
192 rel.get('link', 'unknown')
193 ])
195 print(f"[OK] Detailed files saved for {len(results)} pathways")
198def print_statistics(results):
199 """Print analysis statistics."""
200 print(f"\n{'='*70}")
201 print("PATHWAY ANALYSIS STATISTICS")
202 print(f"{'='*70}")
204 # Total stats
205 total_pathways = len(results)
206 total_interactions = sum(r['num_relations'] for r in results)
207 total_genes = sum(r['num_entries'] for r in results)
209 print(f"\nOverall:")
210 print(f" Total pathways: {total_pathways}")
211 print(f" Total KGML entries (including groups/compounds/maps): {total_genes}")
212 print(f" Total relation subtype records: {total_interactions}")
214 # Largest pathways
215 print(f"\nLargest pathways (by entry count):")
216 sorted_by_size = sorted(results, key=lambda x: x['num_entries'], reverse=True)
217 for i, result in enumerate(sorted_by_size[:10], 1):
218 print(f" {i}. {result['pathway_id']}: {result['num_entries']} entries")
219 print(f" {result['pathway_name']}")
221 # Most connected pathways
222 print(f"\nMost connected pathways (by interactions):")
223 sorted_by_connections = sorted(results, key=lambda x: x['num_relations'], reverse=True)
224 for i, result in enumerate(sorted_by_connections[:10], 1):
225 print(f" {i}. {result['pathway_id']}: {result['num_relations']} interactions")
226 print(f" {result['pathway_name']}")
228 # Interaction type distribution
229 print(f"\nInteraction type distribution:")
230 all_types = Counter()
231 for result in results:
232 for rel_type, count in result['relation_types'].items():
233 all_types[rel_type] += count
235 for rel_type, count in all_types.most_common():
236 percentage = (count / total_interactions) * 100 if total_interactions > 0 else 0
237 print(f" {rel_type}: {count} ({percentage:.1f}%)")
240def main():
241 """Main analysis workflow."""
242 parser = argparse.ArgumentParser(
243 description="Analyze KEGG pathways for an organism",
244 formatter_class=argparse.RawDescriptionHelpFormatter,
245 epilog="""
246Examples:
247 python pathway_analysis.py hsa ./human_pathways
248 python pathway_analysis.py mmu ./mouse_pathways --limit 50
250Organism codes:
251 hsa = Homo sapiens (human)
252 mmu = Mus musculus (mouse)
253 dme = Drosophila melanogaster
254 sce = Saccharomyces cerevisiae (yeast)
255 eco = Escherichia coli
256 """
257 )
258 parser.add_argument("organism", help="KEGG organism code (e.g., hsa, mmu)")
259 parser.add_argument("output_dir", help="Output directory for results")
260 parser.add_argument("--limit", type=int, default=None,
261 help="Limit analysis to first N pathways")
263 args = parser.parse_args()
265 print("=" * 70)
266 print("BIOSERVICES: KEGG Pathway Network Analysis")
267 print("=" * 70)
269 # Create output directory
270 os.makedirs(args.output_dir, exist_ok=True)
272 # Initialize KEGG
273 kegg = KEGG()
274 kegg.services.url = "https://rest.kegg.jp"
276 # Get all pathways
277 pathway_ids = get_all_pathways(kegg, args.organism)
279 if not pathway_ids:
280 print(f"\n[FAIL] No pathways found for {args.organism}")
281 sys.exit(1)
283 # Analyze pathways
284 results = analyze_all_pathways(kegg, pathway_ids, args.limit)
286 if not results:
287 print("\n[FAIL] No pathways successfully analyzed")
288 sys.exit(1)
290 # Print statistics
291 print_statistics(results)
293 # Save results
294 summary_file = os.path.join(args.output_dir, "pathway_summary.csv")
295 save_pathway_summary(results, summary_file)
297 sif_file = os.path.join(args.output_dir, "all_interactions.sif")
298 save_interactions_sif(results, sif_file)
300 save_detailed_pathway_info(results, args.output_dir)
302 # Final summary
303 print(f"\n{'='*70}")
304 print("OUTPUT FILES")
305 print(f"{'='*70}")
306 print(f" Summary: {summary_file}")
307 print(f" Interactions: {sif_file}")
308 print(f" Detailed: {args.output_dir}/pathways/")
309 print(f"{'='*70}")
312if __name__ == "__main__":
313 main()