Skills › Scientific Agent Skills › gget › gene_analysis.py
gene_analysis.py
200 lines · skills/gget/scripts/gene_analysis.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Gene Analysis Script
4Quick analysis of a gene: search, info, sequences, expression, and enrichment
5"""
7import argparse
8import sys
9import gget
12def fasta_text(sequences):
13 """Render what gget.seq returned as FASTA text.
15 gget.seq returns a list of FASTA lines (header, sequence, ...); older
16 releases returned one already-joined string. Writing the list straight to a
17 file raises TypeError, so normalise both shapes here.
18 """
19 if sequences is None:
20 return ""
21 if isinstance(sequences, str):
22 return sequences if sequences.endswith("\n") else sequences + "\n"
23 return "\n".join(sequences) + "\n"
26def analyze_gene(gene_name, species="homo_sapiens", output_prefix=None):
27 """
28 Perform comprehensive analysis of a gene.
30 Args:
31 gene_name: Gene symbol to analyze
32 species: Species name (default: homo_sapiens)
33 output_prefix: Prefix for output files (default: gene_name)
34 """
35 if output_prefix is None:
36 output_prefix = gene_name.lower()
38 print(f"Analyzing gene: {gene_name}")
39 print("=" * 60)
41 # Step 1: Search for the gene
42 print("\n1. Searching for gene...")
43 # search is a substring/description search, not an exact symbol resolver.
44 search_results = gget.search([gene_name], species=species)
46 if search_results is None or search_results.empty:
47 print(f"Error: Gene '{gene_name}' not found in {species}")
48 return False
50 search_results = search_results.loc[
51 search_results["gene_name"].eq(gene_name)
52 ].drop_duplicates("ensembl_id")
53 if len(search_results) != 1:
54 print(f"Error: Expected one exact symbol match for '{gene_name}'; found {len(search_results)}")
55 return False
57 gene_id = search_results["ensembl_id"].iloc[0]
58 print(f" Found: {gene_id}")
59 print(f" Description: {search_results['ensembl_description'].iloc[0]}")
61 # Step 2: Get detailed information
62 print("\n2. Getting detailed information...")
63 gene_info = gget.info([gene_id], pdb=True)
64 if gene_info is None or gene_info.empty:
65 print(f"Error: No annotation returned for {gene_id}")
66 return False
67 gene_info.to_csv(f"{output_prefix}_info.csv", index=True, index_label="query_ensembl_id")
68 print(f" Saved to: {output_prefix}_info.csv")
70 if "uniprot_id" in gene_info.columns and gene_info["uniprot_id"].iloc[0]:
71 print(f" UniProt ID: {gene_info['uniprot_id'].iloc[0]}")
72 if "pdb_id" in gene_info.columns and gene_info["pdb_id"].iloc[0]:
73 print(f" PDB IDs: {gene_info['pdb_id'].iloc[0]}")
75 # Step 3: Get sequences
76 print("\n3. Retrieving sequences...")
77 nucleotide_seq = gget.seq([gene_id])
78 protein_seq = gget.seq([gene_id], translate=True)
79 if not nucleotide_seq or not protein_seq:
80 print("Error: Required sequence retrieval failed; no empty FASTA files written")
81 return False
83 with open(f"{output_prefix}_nucleotide.fasta", "w") as f:
84 f.write(fasta_text(nucleotide_seq))
85 print(f" Nucleotide sequence saved to: {output_prefix}_nucleotide.fasta")
87 with open(f"{output_prefix}_protein.fasta", "w") as f:
88 f.write(fasta_text(protein_seq))
89 print(f" Protein sequence saved to: {output_prefix}_protein.fasta")
91 archs_species = {"human": "human", "homo_sapiens": "human",
92 "mouse": "mouse", "mus_musculus": "mouse"}.get(species)
94 # Step 4: Get tissue expression
95 print("\n4. Getting tissue expression...")
96 try:
97 if archs_species is None:
98 raise ValueError("ARCHS4 tissue mode supports only human and mouse")
99 tissue_expr = gget.archs4(gene_name, which="tissue", species=archs_species)
100 tissue_expr.to_csv(f"{output_prefix}_tissue_expression.csv", index=False)
101 print(f" Saved to: {output_prefix}_tissue_expression.csv")
103 # Show top tissues
104 top_tissues = tissue_expr.nlargest(5, "median")
105 print("\n Top expressing tissues:")
106 for _, row in top_tissues.iterrows():
107 print(f" {row['id']}: median = {row['median']:.2f}")
108 except Exception as e:
109 print(f" Warning: Could not retrieve ARCHS4 data: {e}")
111 # Step 5: Find correlated genes
112 print("\n5. Finding correlated genes...")
113 try:
114 if archs_species != "human":
115 raise ValueError("ARCHS4 correlation mode uses human data")
116 correlated = gget.archs4(gene_name, which="correlation")
117 correlated.to_csv(f"{output_prefix}_correlated_genes.csv", index=False)
118 print(f" Saved to: {output_prefix}_correlated_genes.csv")
120 # Show top correlated
121 print("\n Top 10 correlated genes:")
122 for _, row in correlated.head(10).iterrows():
123 print(f" {row['gene_symbol']}: r = {row['pearson_correlation']:.3f}")
124 except Exception as e:
125 print(f" Warning: Could not retrieve correlation data: {e}")
127 # Step 6: Get disease associations
128 print("\n6. Getting disease associations...")
129 try:
130 if archs_species != "human":
131 raise ValueError("Open Targets requires a human Ensembl gene ID")
132 diseases = gget.opentargets(gene_id, resource="diseases", limit=10)
133 diseases.to_csv(f"{output_prefix}_diseases.csv", index=False)
134 print(f" Saved to: {output_prefix}_diseases.csv")
136 print("\n Top 5 disease associations:")
137 for _, row in diseases.head(5).iterrows():
138 print(f" {row['disease.name']}: score = {row['score']:.3f}")
139 except Exception as e:
140 print(f" Warning: Could not retrieve disease data: {e}")
142 # Step 7: Get drug associations
143 print("\n7. Getting drug associations...")
144 try:
145 if archs_species != "human":
146 raise ValueError("Open Targets requires a human Ensembl gene ID")
147 drugs = gget.opentargets(gene_id, resource="drugs", limit=10)
148 if len(drugs) > 0:
149 drugs.to_csv(f"{output_prefix}_drugs.csv", index=False)
150 print(f" Saved to: {output_prefix}_drugs.csv")
151 print(f"\n Found {len(drugs)} drug associations")
152 else:
153 print(" No drug associations found")
154 except Exception as e:
155 print(f" Warning: Could not retrieve drug data: {e}")
157 print("\n" + "=" * 60)
158 print("Analysis complete!")
159 print(f"\nOutput files (prefix: {output_prefix}):")
160 print(f" - {output_prefix}_info.csv")
161 print(f" - {output_prefix}_nucleotide.fasta")
162 print(f" - {output_prefix}_protein.fasta")
163 print(f" - {output_prefix}_tissue_expression.csv")
164 print(f" - {output_prefix}_correlated_genes.csv")
165 print(f" - {output_prefix}_diseases.csv")
166 print(f" - {output_prefix}_drugs.csv (if available)")
168 return True
171def main():
172 parser = argparse.ArgumentParser(
173 description="Perform comprehensive analysis of a gene using gget"
174 )
175 parser.add_argument("gene", help="Gene symbol to analyze")
176 parser.add_argument(
177 "-s",
178 "--species",
179 default="homo_sapiens",
180 help="Species (default: homo_sapiens)",
181 )
182 parser.add_argument(
183 "-o", "--output", help="Output prefix for files (default: gene name)"
184 )
186 args = parser.parse_args()
188 try:
189 success = analyze_gene(args.gene, args.species, args.output)
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 sys.exit(1)
199if __name__ == "__main__":
200 main()