Skills › Scientific Agent Skills › gget › batch_sequence_analysis.py
batch_sequence_analysis.py
203 lines · skills/gget/scripts/batch_sequence_analysis.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Batch Sequence Analysis Script
4Analyze multiple sequences: BLAST, alignment, and structure prediction
5"""
7import argparse
8import re
9import sys
10from pathlib import Path
11import gget
14def read_fasta(fasta_file):
15 """Read sequences from FASTA file."""
16 sequences = []
17 current_id = None
18 current_seq = []
20 with open(fasta_file, "r") as f:
21 for line in f:
22 line = line.strip()
23 if line.startswith(">"):
24 if current_id:
25 sequences.append({"id": current_id, "seq": "".join(current_seq)})
26 current_id = line[1:]
27 current_seq = []
28 else:
29 current_seq.append(line)
31 if current_id:
32 sequences.append({"id": current_id, "seq": "".join(current_seq)})
34 return sequences
37def analyze_sequences(
38 fasta_file,
39 blast_db="default",
40 align=True,
41 predict_structure=False,
42 output_dir="output",
43):
44 """
45 Perform batch sequence analysis.
47 Args:
48 fasta_file: Path to FASTA file with sequences
49 blast_db: BLAST database to search (default: infer nt/nr from sequence)
50 align: Whether to perform multiple sequence alignment
51 predict_structure: Show the legacy placeholder; no prediction is executed
52 output_dir: Output directory for results
53 """
54 output_path = Path(output_dir)
55 output_path.mkdir(parents=True, exist_ok=True)
57 print(f"Batch Sequence Analysis")
58 print("=" * 60)
59 print(f"Input file: {fasta_file}")
60 print(f"Output directory: {output_dir}")
61 print("")
63 # Read sequences
64 print("Reading sequences...")
65 sequences = read_fasta(fasta_file)
66 print(f"Found {len(sequences)} sequences\n")
68 # Step 1: BLAST each sequence
69 print("Step 1: Running BLAST searches...")
70 print("-" * 60)
71 used_names = set()
72 for i, seq_data in enumerate(sequences):
73 print(f"\n{i+1}. BLASTing {seq_data['id']}...")
74 try:
75 blast_results = gget.blast(
76 seq_data["seq"], database=blast_db, limit=10, save=False
77 )
79 # FASTA headers are data, not file paths. Keep outputs contained and
80 # avoid overwriting repeated or similarly sanitized identifiers.
81 stem = re.sub(r"[^A-Za-z0-9_.-]", "_", seq_data["id"]).strip(".")[:100] or "sequence"
82 name = stem
83 suffix = 2
84 while name in used_names:
85 name = f"{stem}_{suffix}"
86 suffix += 1
87 used_names.add(name)
88 output_file = output_path / f"{name}_blast.csv"
89 blast_results.to_csv(output_file, index=False)
90 print(f" Results saved to: {output_file}")
92 if len(blast_results) > 0:
93 print(f" Top hit: {blast_results.iloc[0]['Description']}")
94 print(
95 f" Max Score: {blast_results.iloc[0]['Max Score']}, "
96 f"Query Coverage: {blast_results.iloc[0]['Query Coverage']}"
97 )
98 except Exception as e:
99 print(f" Error: {e}")
101 # Step 2: Multiple sequence alignment
102 if align and len(sequences) > 1:
103 print("\n\nStep 2: Multiple sequence alignment...")
104 print("-" * 60)
105 try:
106 # gget.muscle returns None: it writes the alignment to `out` (or
107 # prints it when `out` is omitted), so the file must be requested
108 # rather than written from the return value.
109 alignment_file = output_path / "alignment.afa"
110 gget.muscle(fasta_file, out=str(alignment_file))
111 print(f"Alignment saved to: {alignment_file}")
112 except Exception as e:
113 print(f"Error in alignment: {e}")
114 else:
115 print("\n\nStep 2: Skipping alignment (only 1 sequence or disabled)")
117 # Step 3: Structure prediction (optional)
118 if predict_structure:
119 print("\n\nStep 3: Structure prediction placeholder (not executed)")
120 print("-" * 60)
121 print(
122 "The upstream gget AlphaFold wrapper is deprecated and no longer maintained."
123 )
125 for i, seq_data in enumerate(sequences):
126 print(f"\n{i+1}. Predicting structure for {seq_data['id']}...")
127 try:
128 structure_dir = output_path / f"structure_{seq_data['id']}"
129 # Uncomment to run AlphaFold prediction:
130 # gget.alphafold(seq_data['seq'], out=str(structure_dir))
131 # print(f" Structure saved to: {structure_dir}")
132 print(
133 " (Prediction skipped - uncomment code to run AlphaFold prediction)"
134 )
135 except Exception as e:
136 print(f" Error: {e}")
137 else:
138 print("\n\nStep 3: Structure prediction disabled")
140 # Summary
141 print("\n" + "=" * 60)
142 print("Batch analysis complete!")
143 print(f"\nResults saved to: {output_dir}/")
144 print(f" - BLAST results: *_blast.csv")
145 if align and len(sequences) > 1:
146 print(f" - Alignment: alignment.afa")
147 if predict_structure:
148 print(" - Structure prediction was not run; no structure files created")
150 return True
153def main():
154 parser = argparse.ArgumentParser(
155 description="Perform batch sequence analysis using gget"
156 )
157 parser.add_argument("fasta", help="Input FASTA file with sequences")
158 parser.add_argument(
159 "-db",
160 "--database",
161 default="default",
162 help="BLAST database (default: infer nr for proteins, nt for nucleotides)",
163 )
164 parser.add_argument(
165 "--no-align", action="store_true", help="Skip multiple sequence alignment"
166 )
167 parser.add_argument(
168 "--predict-structure",
169 action="store_true",
170 help="Show the legacy structure-prediction placeholder; no prediction is executed",
171 )
172 parser.add_argument(
173 "-o", "--output", default="output", help="Output directory (default: output)"
174 )
176 args = parser.parse_args()
178 if not Path(args.fasta).exists():
179 print(f"Error: File not found: {args.fasta}")
180 sys.exit(1)
182 try:
183 success = analyze_sequences(
184 args.fasta,
185 blast_db=args.database,
186 align=not args.no_align,
187 predict_structure=args.predict_structure,
188 output_dir=args.output,
189 )
190 sys.exit(0 if success else 1)
191 except KeyboardInterrupt:
192 print("\n\nAnalysis interrupted by user")
193 sys.exit(1)
194 except Exception as e:
195 print(f"\n\nError: {e}")
196 import traceback
198 traceback.print_exc()
199 sys.exit(1)
202if __name__ == "__main__":
203 main()