Skills › Scientific Agent Skills › arboreto › basic_grn_inference.py
basic_grn_inference.py
104 lines · skills/arboreto/scripts/basic_grn_inference.py at 5f2bc03
1#!/usr/bin/env python3
2"""Infer candidate TF-target associations from a observations-by-genes TSV.
4Use the exact compatibility environment in SKILL.md. The output is headerless
5TF, target, importance TSV; scores do not establish causal regulation.
6"""
8import argparse
9import csv
10from importlib.metadata import version
12import numpy as np
13import pandas as pd
16def load_expression(expression_file, index_col=None):
17 """Read numeric data without silently renaming duplicate gene identifiers."""
18 with open(expression_file, encoding="utf-8-sig", newline="") as handle:
19 header = next(csv.reader(handle, delimiter="\t"), [])
20 if index_col is not None:
21 if not 0 <= index_col < len(header):
22 raise ValueError("index_col must identify a column in the input file")
23 header = header[:index_col] + header[index_col + 1:]
24 if not header or any(not name.strip() for name in header):
25 raise ValueError("Gene names must be nonempty; use --index-col 0 for row IDs")
26 if len(header) != len(set(header)):
27 raise ValueError("Duplicate gene names: resolve identifiers before inference")
28 frame = pd.read_csv(expression_file, sep="\t", index_col=index_col)
29 if frame.shape[0] < 2 or frame.shape[1] < 2:
30 raise ValueError("At least two observations and two genes are required")
31 try:
32 frame = frame.apply(pd.to_numeric, errors="raise")
33 except (TypeError, ValueError) as exc:
34 raise ValueError("Expression must be numeric; use --index-col 0 for row IDs") from exc
35 if not np.isfinite(frame.to_numpy()).all():
36 raise ValueError("Expression contains missing or nonfinite values")
37 return frame
40def infer_network(expression_data, tf_names, seed, limit, workers):
41 """Use a fresh bounded client; keep legacy Dask configuration explicit."""
42 import dask
44 # PyPI Arboreto 0.1.6 creates an empty metadata graph. This tested Dask
45 # release still offers the legacy backend that accepts that graph.
46 if version("dask") == "2024.7.1":
47 dask.config.set({"dataframe.query-planning": False})
48 from arboreto.algo import grnboost2
49 from distributed import Client, LocalCluster
51 with LocalCluster(
52 n_workers=workers, threads_per_worker=1, dashboard_address=None
53 ) as cluster, Client(cluster) as client:
54 return grnboost2(
55 expression_data=expression_data, tf_names=tf_names, seed=seed,
56 limit=limit, client_or_address=client, verbose=True,
57 )
60def run_grn_inference(expression_file, output_file, tf_file=None, seed=777,
61 limit=None, index_col=None, workers=1):
62 """Validate inputs, infer candidates, and write a headerless adjacency TSV."""
63 if limit is not None and limit <= 0:
64 raise ValueError("limit must be a positive integer")
65 if workers <= 0:
66 raise ValueError("workers must be a positive integer")
67 expression_data = load_expression(expression_file, index_col=index_col)
68 print(f"[OK] Loaded {expression_data.shape[0]} observations x {expression_data.shape[1]} genes")
69 tf_names = "all"
70 if tf_file:
71 with open(tf_file, encoding="utf-8") as handle:
72 requested = list(dict.fromkeys(line.strip() for line in handle if line.strip()))
73 tf_names = [name for name in requested if name in expression_data.columns]
74 if not tf_names:
75 raise ValueError("TF file has no names matching the expression genes")
76 print(f"[OK] Matched {len(tf_names)} of {len(requested)} TF names")
77 else:
78 print("[OK] All genes are candidate regulators; the TF column is not restricted to known TFs")
79 network = infer_network(expression_data, tf_names, seed, limit, workers)
80 if network.empty:
81 raise RuntimeError("Inference returned no links; inspect worker warnings and input variation")
82 if not np.isfinite(network["importance"]).all():
83 raise RuntimeError("Inference returned nonfinite importance scores")
84 network.to_csv(output_file, sep="\t", index=False, header=False)
85 print(f"[OK] Saved {len(network)} candidate links to {output_file}")
86 print(network.head(10).to_string(index=False))
87 return network
90def main():
91 parser = argparse.ArgumentParser(description=__doc__)
92 parser.add_argument("expression_file", help="Numeric TSV: observations in rows, genes in columns")
93 parser.add_argument("output_file", help="Headerless TF, target, importance TSV")
94 parser.add_argument("--tf-file", help="TF names, one per line")
95 parser.add_argument("--seed", type=int, default=777, help="Regressor seed (default: 777)")
96 parser.add_argument("--limit", type=int, help="Return the top N links globally; positive integer")
97 parser.add_argument("--index-col", type=int, help="Zero-based row identifier column, typically 0")
98 parser.add_argument("--workers", type=int, default=1, help="Local worker processes (default: 1)")
99 args = parser.parse_args()
100 run_grn_inference(**vars(args))
103if __name__ == "__main__":
104 main()