Skills › Scientific Agent Skills › pysam › inspect_hts.py
inspect_hts.py
497 lines · skills/pysam/scripts/inspect_hts.py at 5f2bc03
1#!/usr/bin/env python3
2"""Inspect a local HTS/sequence file and emit a bounded JSON summary.
4The default report reads headers and index metadata only. It does not emit
5alignment query names, FASTX record names, VCF sample names, or full headers.
6CRAM requires an explicit local --reference to avoid hidden reference lookup.
7"""
9from __future__ import annotations
11import argparse
12import json
13import sys
14from pathlib import Path
15from typing import Any, Optional
17import pysam
20KINDS = ("alignment", "variant", "fasta", "fastx", "tabix")
23def positive_int(value: str) -> int:
24 parsed = int(value)
25 if parsed <= 0:
26 raise argparse.ArgumentTypeError("must be greater than zero")
27 return parsed
30def nonnegative_int(value: str) -> int:
31 parsed = int(value)
32 if parsed < 0:
33 raise argparse.ArgumentTypeError("must be zero or greater")
34 return parsed
37def detect_kind(path: Path) -> str:
38 name = path.name.lower()
40 if name.endswith((".bam", ".sam", ".cram")):
41 return "alignment"
42 if name.endswith((".vcf", ".vcf.gz", ".vcf.bgz", ".bcf")):
43 return "variant"
44 if name.endswith(
45 (
46 ".fa",
47 ".fasta",
48 ".fna",
49 ".fas",
50 ".fa.gz",
51 ".fasta.gz",
52 ".fna.gz",
53 ".fas.gz",
54 )
55 ):
56 return "fasta"
57 if name.endswith(
58 (
59 ".fq",
60 ".fastq",
61 ".fq.gz",
62 ".fastq.gz",
63 ".fq.bgz",
64 ".fastq.bgz",
65 )
66 ):
67 return "fastx"
68 if name.endswith(
69 (
70 ".bed.gz",
71 ".bed.bgz",
72 ".gff.gz",
73 ".gff3.gz",
74 ".gtf.gz",
75 ".tsv.gz",
76 ".tab.gz",
77 )
78 ):
79 return "tabix"
81 raise ValueError(
82 "could not infer file kind from suffix; pass --kind "
83 + ", ".join(KINDS)
84 )
87def build_parser() -> argparse.ArgumentParser:
88 parser = argparse.ArgumentParser(
89 description=__doc__,
90 formatter_class=argparse.ArgumentDefaultsHelpFormatter,
91 )
92 parser.add_argument("input", type=Path, help="local genomic data file")
93 parser.add_argument(
94 "--kind",
95 choices=("auto",) + KINDS,
96 default="auto",
97 help="file interface to use",
98 )
99 parser.add_argument(
100 "--reference",
101 type=Path,
102 help="local indexed reference FASTA; required for CRAM",
103 )
104 parser.add_argument(
105 "--index",
106 type=Path,
107 help="explicit local BAI/CSI/CRAI/TBI/FAI index",
108 )
109 parser.add_argument(
110 "--threads",
111 type=positive_int,
112 default=1,
113 help="HTSlib decompression threads",
114 )
115 parser.add_argument(
116 "--max-items",
117 type=nonnegative_int,
118 default=100,
119 help="maximum contigs and metadata keys to include; 0 includes all",
120 )
121 parser.add_argument(
122 "--include-identifiers",
123 action="store_true",
124 help="include VCF sample names and alignment read-group IDs",
125 )
126 parser.add_argument(
127 "--output",
128 type=Path,
129 help="write JSON to a new file instead of stdout; refuses overwrite",
130 )
131 return parser
134def require_local_file(path: Path, label: str) -> Path:
135 resolved = path.expanduser()
136 if not resolved.exists():
137 raise FileNotFoundError(f"{label} does not exist: {resolved}")
138 if not resolved.is_file():
139 raise ValueError(f"{label} is not a regular file: {resolved}")
140 return resolved
143def limited(values: list[Any], maximum: int) -> tuple[list[Any], bool]:
144 if maximum == 0 or len(values) <= maximum:
145 return values, False
146 return values[:maximum], True
149def contig_records(
150 names: tuple[str, ...] | list[str],
151 lengths: tuple[int, ...] | list[int],
152 maximum: int,
153) -> tuple[list[dict[str, Any]], bool]:
154 records = [
155 {"name": name, "length": int(length)}
156 for name, length in zip(names, lengths)
157 ]
158 return limited(records, maximum)
161def alignment_mode(path: Path) -> str:
162 name = path.name.lower()
163 if name.endswith(".sam"):
164 return "r"
165 if name.endswith(".cram"):
166 return "rc"
167 return "rb"
170def inspect_alignment(
171 path: Path,
172 *,
173 reference: Optional[Path],
174 index: Optional[Path],
175 threads: int,
176 max_items: int,
177 include_identifiers: bool,
178) -> dict[str, Any]:
179 if path.name.lower().endswith(".cram") and reference is None:
180 raise ValueError("CRAM inspection requires --reference")
182 kwargs: dict[str, Any] = {"threads": threads}
183 if reference is not None:
184 kwargs["reference_filename"] = str(reference)
185 if index is not None:
186 kwargs["index_filename"] = str(index)
188 with pysam.AlignmentFile(
189 str(path),
190 alignment_mode(path),
191 **kwargs,
192 ) as alignments:
193 contigs, contigs_truncated = contig_records(
194 list(alignments.references),
195 list(alignments.lengths),
196 max_items,
197 )
198 header = alignments.header.to_dict()
199 read_groups = header.get("RG", [])
200 programs = header.get("PG", [])
202 try:
203 has_index = bool(alignments.has_index())
204 except (AttributeError, OSError, ValueError):
205 has_index = False
207 report: dict[str, Any] = {
208 "format": (
209 "CRAM"
210 if alignments.is_cram
211 else "BAM"
212 if alignments.is_bam
213 else "SAM"
214 ),
215 "is_remote": bool(alignments.is_remote),
216 "has_index": has_index,
217 "sort_order": header.get("HD", {}).get("SO"),
218 "reference_count": len(alignments.references),
219 "contigs": contigs,
220 "contigs_truncated": contigs_truncated,
221 "read_group_count": len(read_groups),
222 "program_record_count": len(programs),
223 "comment_count": len(header.get("CO", [])),
224 }
226 if include_identifiers:
227 read_group_ids = [
228 group.get("ID")
229 for group in read_groups
230 if group.get("ID") is not None
231 ]
232 report["read_group_ids"], report["read_group_ids_truncated"] = (
233 limited(read_group_ids, max_items)
234 )
236 if has_index and alignments.is_cram:
237 report["index_statistics"] = None
238 report["index_statistics_note"] = (
239 "CRAI has no mapped/unmapped counts; use a sequential QC scan."
240 )
241 elif has_index:
242 try:
243 statistics = alignments.get_index_statistics()
244 report["index_statistics"] = {
245 "mapped": sum(item.mapped for item in statistics),
246 "unmapped_with_coordinates": sum(
247 item.unmapped for item in statistics
248 ),
249 "no_coordinate": int(alignments.nocoordinate),
250 }
251 except (AttributeError, OSError, ValueError):
252 report["index_statistics"] = None
254 return report
257def variant_index_loaded(variants: pysam.VariantFile) -> bool:
258 try:
259 return variants.index is not None
260 except (AttributeError, OSError, ValueError):
261 return False
264def inspect_variant(
265 path: Path,
266 *,
267 index: Optional[Path],
268 threads: int,
269 max_items: int,
270 include_identifiers: bool,
271) -> dict[str, Any]:
272 kwargs: dict[str, Any] = {"threads": threads}
273 if index is not None:
274 kwargs["index_filename"] = str(index)
276 with pysam.VariantFile(str(path), **kwargs) as variants:
277 header = variants.header
278 contig_values = [
279 {"name": name, "length": metadata.length}
280 for name, metadata in header.contigs.items()
281 ]
282 contigs, contigs_truncated = limited(contig_values, max_items)
283 info_keys, info_truncated = limited(list(header.info), max_items)
284 format_keys, formats_truncated = limited(
285 list(header.formats), max_items
286 )
287 filter_keys, filters_truncated = limited(
288 list(header.filters), max_items
289 )
291 report: dict[str, Any] = {
292 "format": (
293 "BCF"
294 if bool(getattr(variants, "is_bcf", False))
295 else "VCF"
296 ),
297 "is_remote": bool(variants.is_remote),
298 "has_index": variant_index_loaded(variants),
299 "contig_count": len(header.contigs),
300 "contigs": contigs,
301 "contigs_truncated": contigs_truncated,
302 "sample_count": len(header.samples),
303 "info_keys": info_keys,
304 "info_keys_truncated": info_truncated,
305 "format_keys": format_keys,
306 "format_keys_truncated": formats_truncated,
307 "filter_keys": filter_keys,
308 "filter_keys_truncated": filters_truncated,
309 }
310 if include_identifiers:
311 sample_names, samples_truncated = limited(
312 list(header.samples), max_items
313 )
314 report["sample_names"] = sample_names
315 report["sample_names_truncated"] = samples_truncated
316 return report
319def inspect_fasta(
320 path: Path,
321 *,
322 index: Optional[Path],
323 max_items: int,
324) -> dict[str, Any]:
325 # An implicit FastaFile index can be created on open. An explicit index
326 # disables that write and preserves this inspector's read-only contract.
327 fasta_index = require_local_file(
328 index if index is not None else Path(str(path) + ".fai"),
329 "FASTA index (create explicitly with pysam.faidx first)",
330 )
331 kwargs: dict[str, Any] = {"filepath_index": str(fasta_index)}
333 with pysam.FastaFile(str(path), **kwargs) as fasta:
334 contigs, contigs_truncated = contig_records(
335 list(fasta.references),
336 list(fasta.lengths),
337 max_items,
338 )
339 return {
340 "format": "FASTA",
341 "random_access": True,
342 "reference_count": fasta.nreferences,
343 "total_bases": sum(fasta.lengths),
344 "contigs": contigs,
345 "contigs_truncated": contigs_truncated,
346 }
349def inspect_fastx(path: Path) -> dict[str, Any]:
350 with pysam.FastxFile(str(path), persist=False):
351 return {
352 "format": "FASTA/FASTQ",
353 "random_access": False,
354 "content_scanned": False,
355 "note": (
356 "FastxFile has no header index; use a streaming QC workflow "
357 "to summarize records."
358 ),
359 }
362def inspect_tabix(
363 path: Path,
364 *,
365 index: Optional[Path],
366 threads: int,
367 max_items: int,
368) -> dict[str, Any]:
369 # TabixFile defaults to .tbi only; CSI requires an explicit index path.
370 kwargs: dict[str, Any] = {"threads": threads}
371 if index is not None:
372 kwargs["index"] = str(index)
373 elif not Path(str(path) + ".tbi").exists() and Path(str(path) + ".csi").is_file():
374 kwargs["index"] = str(path) + ".csi"
376 with pysam.TabixFile(str(path), **kwargs) as table:
377 contigs, contigs_truncated = limited(
378 list(table.contigs), max_items
379 )
380 header_line_count = sum(1 for _ in table.header)
381 return {
382 "format": "TABIX",
383 "has_index": True,
384 "contig_count": len(table.contigs),
385 "contigs": contigs,
386 "contigs_truncated": contigs_truncated,
387 "header_line_count": header_line_count,
388 }
391def inspect_file(args: argparse.Namespace) -> dict[str, Any]:
392 path = require_local_file(args.input, "input")
393 reference = (
394 require_local_file(args.reference, "reference")
395 if args.reference is not None
396 else None
397 )
398 index = (
399 require_local_file(args.index, "index")
400 if args.index is not None
401 else None
402 )
403 kind = detect_kind(path) if args.kind == "auto" else args.kind
405 if kind == "alignment":
406 details = inspect_alignment(
407 path,
408 reference=reference,
409 index=index,
410 threads=args.threads,
411 max_items=args.max_items,
412 include_identifiers=args.include_identifiers,
413 )
414 elif kind == "variant":
415 details = inspect_variant(
416 path,
417 index=index,
418 threads=args.threads,
419 max_items=args.max_items,
420 include_identifiers=args.include_identifiers,
421 )
422 elif kind == "fasta":
423 details = inspect_fasta(
424 path,
425 index=index,
426 max_items=args.max_items,
427 )
428 elif kind == "fastx":
429 details = inspect_fastx(path)
430 else:
431 details = inspect_tabix(
432 path,
433 index=index,
434 threads=args.threads,
435 max_items=args.max_items,
436 )
438 return {
439 "schema_version": "1.0",
440 "pysam_version": pysam.__version__,
441 "samtools_version": getattr(
442 pysam, "__samtools_version__", None
443 ),
444 "file_name": path.name,
445 "size_bytes": path.stat().st_size,
446 "kind": kind,
447 "identifiers_included": args.include_identifiers,
448 "details": details,
449 }
452def write_report(report: dict[str, Any], output: Optional[Path]) -> None:
453 text = json.dumps(
454 report,
455 indent=2,
456 ensure_ascii=False,
457 allow_nan=False,
458 )
459 if output is None:
460 sys.stdout.write(text + "\n")
461 return
463 destination = output.expanduser()
464 with destination.open("x", encoding="utf-8") as handle:
465 handle.write(text + "\n")
468def main(argv: Optional[list[str]] = None) -> int:
469 parser = build_parser()
470 args = parser.parse_args(argv)
472 if args.output is not None:
473 try:
474 if args.output.expanduser().resolve() == args.input.expanduser().resolve():
475 parser.error("--output must not overwrite the input file")
476 except OSError:
477 pass
479 try:
480 report = inspect_file(args)
481 write_report(report, args.output)
482 except (
483 FileExistsError,
484 FileNotFoundError,
485 IndexError,
486 KeyError,
487 OSError,
488 TypeError,
489 ValueError,
490 ) as error:
491 print(f"inspect_hts: {error}", file=sys.stderr)
492 return 2
493 return 0
496if __name__ == "__main__":
497 raise SystemExit(main())