Skills › Scientific Agent Skills › gtars › coverage_preflight.py
coverage_preflight.py
262 lines · skills/gtars/scripts/coverage_preflight.py at 5f2bc03
1#!/usr/bin/env python3
2"""Preflight local Gtars uniwig coverage and bigWig generation."""
4from __future__ import annotations
6import argparse
7import math
8import sys
10from _common import (
11 HARD_MAX_BYTES,
12 HARD_MAX_RECORDS,
13 HARD_MAX_WORKERS,
14 SafetyError,
15 add_path_mode_argument,
16 display_path,
17 fail_json,
18 inspect_bed,
19 int_type,
20 load_chrom_sizes,
21 local_path,
22 print_json,
23 prefix_collision_count,
24 sha256_file,
25)
28TOOL = "gtars-coverage-preflight"
31def build_parser() -> argparse.ArgumentParser:
32 parser = argparse.ArgumentParser(
33 description=(
34 "Check a sorted local BED/narrowPeak input, chromosome sizes, bounds, "
35 "and output budget for gtars uniwig 0.10.0. Nothing is generated."
36 )
37 )
38 parser.add_argument("--input", required=True)
39 parser.add_argument("--input-type", choices=("bed", "narrowpeak"), default="bed")
40 parser.add_argument("--chrom-sizes", required=True)
41 parser.add_argument("--assembly", required=True)
42 parser.add_argument("--output-prefix", required=True)
43 parser.add_argument(
44 "--output-type",
45 choices=("bw", "wig", "bedgraph", "npy"),
46 default="bw",
47 )
48 parser.add_argument(
49 "--count-type",
50 choices=("start", "end", "core", "all"),
51 default="core",
52 )
53 parser.add_argument(
54 "--smooth-size",
55 type=int_type(minimum=1, maximum=10_000_000, label="smooth-size"),
56 default=5,
57 )
58 parser.add_argument(
59 "--step-size",
60 type=int_type(minimum=1, maximum=10_000_000, label="step-size"),
61 default=1,
62 )
63 parser.add_argument(
64 "--threads",
65 type=int_type(minimum=1, maximum=HARD_MAX_WORKERS, label="threads"),
66 default=1,
67 )
68 parser.add_argument(
69 "--zoom",
70 type=int_type(minimum=0, maximum=10_000, label="zoom"),
71 default=1,
72 )
73 parser.add_argument(
74 "--streaming",
75 action="store_true",
76 help="Use BED streaming; only wig/bedgraph output is supported upstream.",
77 )
78 parser.add_argument(
79 "--dense",
80 type=int_type(minimum=-1, maximum=MAX_DENSE_GAP, label="dense"),
81 default=100,
82 help="Streaming gap fill: 0 sparse, -1 fully dense, N fills gaps <=N.",
83 )
84 parser.add_argument(
85 "--max-bytes",
86 type=int_type(minimum=1, maximum=HARD_MAX_BYTES, label="max-bytes"),
87 default=2 * 1024**3,
88 )
89 parser.add_argument(
90 "--max-records",
91 type=int_type(
92 minimum=1,
93 maximum=HARD_MAX_RECORDS,
94 label="max-records",
95 ),
96 default=5_000_000,
97 )
98 parser.add_argument(
99 "--max-estimated-bytes",
100 type=int_type(
101 minimum=1,
102 maximum=HARD_MAX_BYTES,
103 label="max-estimated-bytes",
104 ),
105 default=4 * 1024**3,
106 help="Planning cap for dense uncompressed coverage values.",
107 )
108 add_path_mode_argument(parser)
109 return parser
112MAX_DENSE_GAP = 10_000_000
115def preflight(args: argparse.Namespace) -> tuple[dict, int]:
116 assembly = args.assembly.strip()
117 if not assembly or len(assembly) > 200:
118 raise SafetyError("assembly must contain 1-200 characters")
119 input_path = local_path(args.input, kind="file")
120 chrom_path = local_path(args.chrom_sizes, kind="file")
121 output_prefix = local_path(args.output_prefix, must_exist=False)
123 chrom_sizes, chrom_order = load_chrom_sizes(
124 chrom_path,
125 max_bytes=min(args.max_bytes, 128 * 1024**2),
126 max_records=min(args.max_records, 1_000_000),
127 )
128 chrom_digest, chrom_bytes = sha256_file(
129 chrom_path,
130 max_bytes=min(args.max_bytes, 128 * 1024**2),
131 )
132 bed = inspect_bed(
133 input_path,
134 max_bytes=args.max_bytes,
135 max_records=args.max_records,
136 chrom_sizes=chrom_sizes,
137 chrom_order=chrom_order,
138 )
139 errors: dict[str, int] = dict(bed["errors"])
140 warnings: dict[str, int] = {}
141 if bed["records"] == 0:
142 errors["no_data_records"] = 1
143 if bed["out_of_order_records"]:
144 errors["sorting_required"] = bed["out_of_order_records"]
145 if args.input_type == "narrowpeak" and (bed["minimum_columns"] or 0) < 10:
146 errors["narrowpeak_requires_ten_columns"] = 1
147 if output_prefix.exists():
148 errors["output_prefix_already_exists"] = 1
149 collisions = prefix_collision_count(output_prefix)
150 if collisions:
151 errors["output_prefix_collisions"] = collisions
152 if args.streaming and args.output_type not in {"wig", "bedgraph"}:
153 errors["streaming_output_type_unsupported"] = 1
154 if args.streaming and args.input_type != "bed":
155 errors["streaming_requires_bed_input"] = 1
157 count_outputs = 3 if args.count_type == "all" else 1
158 genome_span = sum(chrom_sizes.values())
159 estimated_values = sum(
160 math.ceil(length / args.step_size) for length in chrom_sizes.values()
161 ) * count_outputs
162 estimated_uncompressed_bytes = estimated_values * 8
163 if estimated_uncompressed_bytes > args.max_estimated_bytes:
164 errors["estimated_dense_coverage_exceeds_budget"] = 1
165 if args.threads > 1:
166 warnings["parallel_memory_may_scale_with_threads"] = args.threads
167 if args.output_type == "bw":
168 warnings["bigwig_size_is_data_dependent_estimate_is_not_file_size"] = 1
169 if args.streaming and args.dense == -1:
170 warnings["fully_dense_streaming_requested"] = 1
172 argv = [
173 "gtars",
174 "uniwig",
175 "--file",
176 "<input>",
177 "--filetype",
178 args.input_type,
179 "--chromref",
180 "<chrom-sizes>",
181 "--smoothsize",
182 str(args.smooth_size),
183 "--stepsize",
184 str(args.step_size),
185 "--fileheader",
186 "<output-prefix>",
187 "--outputtype",
188 args.output_type,
189 "--counttype",
190 args.count_type,
191 ]
192 if args.streaming:
193 argv.extend(["--streaming", "--dense", str(args.dense)])
194 else:
195 argv.extend(["--threads", str(args.threads), "--zoom", str(args.zoom)])
197 report = {
198 "ok": not errors,
199 "ready_to_execute": not errors,
200 "tool": TOOL,
201 "contract": {
202 "assembly": assembly,
203 "coordinate_system": "0-based-half-open",
204 "gtars_cli_version": "0.10.0",
205 "commands_executed": False,
206 "files_written": False,
207 "network_used": False,
208 "symlinks_allowed": False,
209 },
210 "input": {
211 "path": display_path(input_path, 1, args.path_mode),
212 "sha256": bed["sha256"],
213 "size_bytes": bed["size_bytes"],
214 "records": bed["records"],
215 "minimum_columns": bed["minimum_columns"],
216 "out_of_order_records": bed["out_of_order_records"],
217 },
218 "chromosome_sizes": {
219 "path": display_path(chrom_path, 2, args.path_mode),
220 "sha256": chrom_digest,
221 "size_bytes": chrom_bytes,
222 "contig_count": len(chrom_sizes),
223 "total_span": genome_span,
224 },
225 "output": {
226 "prefix": display_path(output_prefix, 3, args.path_mode),
227 "type": args.output_type,
228 "count_type": args.count_type,
229 },
230 "resource_estimate": {
231 "step_size": args.step_size,
232 "count_outputs": count_outputs,
233 "dense_value_upper_bound": estimated_values,
234 "uncompressed_bytes_proxy": estimated_uncompressed_bytes,
235 "budget_bytes": args.max_estimated_bytes,
236 "not_a_bigwig_size_prediction": True,
237 },
238 "errors": dict(sorted(errors.items())),
239 "warnings": dict(sorted(warnings.items())),
240 "argv_template": argv,
241 "approval_gate": [
242 "review the sorted/bounds report and assembly checksum",
243 "pilot one bounded contig or small synthetic BED first",
244 "set CPU, RAM, disk, and wall-time limits before running",
245 "inspect the resulting bigWig header and chromosome dictionary",
246 ],
247 }
248 return report, 0 if not errors else 2
251def main() -> int:
252 args = build_parser().parse_args()
253 try:
254 report, status = preflight(args)
255 print_json(report)
256 return status
257 except (OSError, SafetyError, UnicodeError) as exc:
258 return fail_json(TOOL, exc)
261if __name__ == "__main__":
262 sys.exit(main())